{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gaussian Process Latent Variable Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [Gaussian Process Latent Variable Model](https://en.wikipedia.org/wiki/Nonlinear_dimensionality_reduction#Gaussian_process_latent_variable_models) (GPLVM) is a dimensionality reduction method that uses a Gaussian process to learn a low-dimensional representation of (potentially) high-dimensional data. In the typical setting of Gaussian process regression, where we are given inputs $X$ and outputs $y$, we choose a kernel and learn hyperparameters that best describe the mapping from $X$ to $y$. In the GPLVM, we are not given $X$: we are only given $y$. So we need to learn $X$ along with the kernel hyperparameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We do not do maximum likelihood inference on $X$. Instead, we set a Gaussian prior for $X$ and learn the mean and variance of the approximate (gaussian) posterior $q(X|y)$. In this notebook, we show how this can be done using the `pyro.contrib.gp` module. In particular we reproduce a result described in [2]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import torch\n",
    "from torch.nn import Parameter\n",
    "\n",
    "import pyro\n",
    "import pyro.contrib.gp as gp\n",
    "import pyro.distributions as dist\n",
    "import pyro.ops.stats as stats\n",
    "\n",
    "smoke_test = ('CI' in os.environ)  # ignore; used to check code integrity in the Pyro repo\n",
    "assert pyro.__version__.startswith('1.9.1')\n",
    "pyro.set_rng_seed(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data we are going to use consists of [single-cell](https://en.wikipedia.org/wiki/Single-cell_analysis) [qPCR](https://en.wikipedia.org/wiki/Real-time_polymerase_chain_reaction) data for 48 genes obtained from mice (Guo *et al.*, [1]). This data is available at the [Open Data Science repository](https://github.com/sods/ods). The data contains 48 columns, with each column corresponding to (normalized) measurements of each gene. Cells differentiate during their development and these data were obtained at various stages of development. The various stages are labelled from the 1-cell stage to the 64-cell stage. For the 32-cell stage, the data is further differentiated into 'trophectoderm' (TE) and 'inner cell mass' (ICM). ICM further differentiates into 'epiblast' (EPI) and 'primitive endoderm' (PE) at the 64-cell stage. Each of the rows in the dataset is labelled with one of these stages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data shape: (437, 48)\n",
      "---------------------\n",
      "\n",
      "Data labels: ['1', '2', '4', '8', '16', '32 TE', '32 ICM', '64 PE', '64 TE', '64 EPI']\n",
      "--------------------------------------------------------------------------------------\n",
      "\n",
      "Show a small subset of the data:\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Actb</th>\n",
       "      <th>Ahcy</th>\n",
       "      <th>Aqp3</th>\n",
       "      <th>Atp12a</th>\n",
       "      <th>Bmp4</th>\n",
       "      <th>Cdx2</th>\n",
       "      <th>Creb312</th>\n",
       "      <th>Cebpa</th>\n",
       "      <th>Dab2</th>\n",
       "      <th>DppaI</th>\n",
       "      <th>...</th>\n",
       "      <th>Sox2</th>\n",
       "      <th>Sall4</th>\n",
       "      <th>Sox17</th>\n",
       "      <th>Snail</th>\n",
       "      <th>Sox13</th>\n",
       "      <th>Tcfap2a</th>\n",
       "      <th>Tcfap2c</th>\n",
       "      <th>Tcf23</th>\n",
       "      <th>Utf1</th>\n",
       "      <th>Tspan8</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.541050</td>\n",
       "      <td>-1.203007</td>\n",
       "      <td>1.030746</td>\n",
       "      <td>1.064808</td>\n",
       "      <td>0.494782</td>\n",
       "      <td>-0.167143</td>\n",
       "      <td>-1.369092</td>\n",
       "      <td>1.083061</td>\n",
       "      <td>0.668057</td>\n",
       "      <td>-1.553758</td>\n",
       "      <td>...</td>\n",
       "      <td>-1.351757</td>\n",
       "      <td>-1.793476</td>\n",
       "      <td>0.783185</td>\n",
       "      <td>-1.408063</td>\n",
       "      <td>-0.031991</td>\n",
       "      <td>-0.351257</td>\n",
       "      <td>-1.078982</td>\n",
       "      <td>0.942981</td>\n",
       "      <td>1.348892</td>\n",
       "      <td>-1.051999</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.680832</td>\n",
       "      <td>-1.355306</td>\n",
       "      <td>2.456375</td>\n",
       "      <td>1.234350</td>\n",
       "      <td>0.645494</td>\n",
       "      <td>1.003868</td>\n",
       "      <td>-1.207595</td>\n",
       "      <td>1.208023</td>\n",
       "      <td>0.800388</td>\n",
       "      <td>-1.435306</td>\n",
       "      <td>...</td>\n",
       "      <td>-1.363533</td>\n",
       "      <td>-1.782172</td>\n",
       "      <td>1.532477</td>\n",
       "      <td>-1.361172</td>\n",
       "      <td>-0.501715</td>\n",
       "      <td>1.082362</td>\n",
       "      <td>-0.930112</td>\n",
       "      <td>1.064399</td>\n",
       "      <td>1.469397</td>\n",
       "      <td>-0.996275</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.056038</td>\n",
       "      <td>-1.280447</td>\n",
       "      <td>2.046133</td>\n",
       "      <td>1.439795</td>\n",
       "      <td>0.828121</td>\n",
       "      <td>0.983404</td>\n",
       "      <td>-1.460032</td>\n",
       "      <td>1.359447</td>\n",
       "      <td>0.530701</td>\n",
       "      <td>-1.340283</td>\n",
       "      <td>...</td>\n",
       "      <td>-1.296802</td>\n",
       "      <td>-1.567402</td>\n",
       "      <td>3.194157</td>\n",
       "      <td>-1.301777</td>\n",
       "      <td>-0.445219</td>\n",
       "      <td>0.031284</td>\n",
       "      <td>-1.005767</td>\n",
       "      <td>1.211529</td>\n",
       "      <td>1.615421</td>\n",
       "      <td>-0.651393</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.732331</td>\n",
       "      <td>-1.326911</td>\n",
       "      <td>2.464234</td>\n",
       "      <td>1.244323</td>\n",
       "      <td>0.654359</td>\n",
       "      <td>0.947023</td>\n",
       "      <td>-1.265609</td>\n",
       "      <td>1.215373</td>\n",
       "      <td>0.765212</td>\n",
       "      <td>-1.431401</td>\n",
       "      <td>...</td>\n",
       "      <td>-1.684100</td>\n",
       "      <td>-1.915556</td>\n",
       "      <td>2.962515</td>\n",
       "      <td>-1.349710</td>\n",
       "      <td>1.875957</td>\n",
       "      <td>1.699892</td>\n",
       "      <td>-1.059458</td>\n",
       "      <td>1.071541</td>\n",
       "      <td>1.476485</td>\n",
       "      <td>-0.699586</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.629333</td>\n",
       "      <td>-1.244308</td>\n",
       "      <td>1.316815</td>\n",
       "      <td>1.304162</td>\n",
       "      <td>0.707552</td>\n",
       "      <td>1.429070</td>\n",
       "      <td>-0.895578</td>\n",
       "      <td>-0.007785</td>\n",
       "      <td>0.644606</td>\n",
       "      <td>-1.381937</td>\n",
       "      <td>...</td>\n",
       "      <td>-1.304653</td>\n",
       "      <td>-1.761825</td>\n",
       "      <td>1.265379</td>\n",
       "      <td>-1.320533</td>\n",
       "      <td>-0.609864</td>\n",
       "      <td>0.413826</td>\n",
       "      <td>-0.888624</td>\n",
       "      <td>1.114394</td>\n",
       "      <td>1.519017</td>\n",
       "      <td>-0.798985</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 48 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "       Actb      Ahcy      Aqp3    Atp12a      Bmp4      Cdx2   Creb312  \\\n",
       "1  0.541050 -1.203007  1.030746  1.064808  0.494782 -0.167143 -1.369092   \n",
       "1  0.680832 -1.355306  2.456375  1.234350  0.645494  1.003868 -1.207595   \n",
       "1  1.056038 -1.280447  2.046133  1.439795  0.828121  0.983404 -1.460032   \n",
       "1  0.732331 -1.326911  2.464234  1.244323  0.654359  0.947023 -1.265609   \n",
       "1  0.629333 -1.244308  1.316815  1.304162  0.707552  1.429070 -0.895578   \n",
       "\n",
       "      Cebpa      Dab2     DppaI  ...      Sox2     Sall4     Sox17     Snail  \\\n",
       "1  1.083061  0.668057 -1.553758  ... -1.351757 -1.793476  0.783185 -1.408063   \n",
       "1  1.208023  0.800388 -1.435306  ... -1.363533 -1.782172  1.532477 -1.361172   \n",
       "1  1.359447  0.530701 -1.340283  ... -1.296802 -1.567402  3.194157 -1.301777   \n",
       "1  1.215373  0.765212 -1.431401  ... -1.684100 -1.915556  2.962515 -1.349710   \n",
       "1 -0.007785  0.644606 -1.381937  ... -1.304653 -1.761825  1.265379 -1.320533   \n",
       "\n",
       "      Sox13   Tcfap2a   Tcfap2c     Tcf23      Utf1    Tspan8  \n",
       "1 -0.031991 -0.351257 -1.078982  0.942981  1.348892 -1.051999  \n",
       "1 -0.501715  1.082362 -0.930112  1.064399  1.469397 -0.996275  \n",
       "1 -0.445219  0.031284 -1.005767  1.211529  1.615421 -0.651393  \n",
       "1  1.875957  1.699892 -1.059458  1.071541  1.476485 -0.699586  \n",
       "1 -0.609864  0.413826 -0.888624  1.114394  1.519017 -0.798985  \n",
       "\n",
       "[5 rows x 48 columns]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# license: Copyright (c) 2014, the Open Data Science Initiative\n",
    "# license: https://www.elsevier.com/legal/elsevier-website-terms-and-conditions\n",
    "URL = \"https://raw.githubusercontent.com/sods/ods/master/datasets/guo_qpcr.csv\"\n",
    "\n",
    "df = pd.read_csv(URL, index_col=0)\n",
    "print(\"Data shape: {}\\n{}\\n\".format(df.shape, \"-\" * 21))\n",
    "print(\"Data labels: {}\\n{}\\n\".format(df.index.unique().tolist(), \"-\" * 86))\n",
    "print(\"Show a small subset of the data:\")\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Modelling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we need to define the output tensor $y$. To predict values for all $48$ genes, we need $48$ Gaussian processes. So the required shape for $y$ is `num_GPs x num_data = 48 x 437`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = torch.tensor(df.values, dtype=torch.get_default_dtype())\n",
    "# we need to transpose data to correct its shape\n",
    "y = data.t()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now comes the most interesting part. We know that the observed data $y$ has latent structure: in particular different datapoints correspond to different cell stages. We would like our GPLVM to learn this structure in an unsupervised manner. In principle, if we do a good job of inference then we should be able to discover this structure---at least if we choose reasonable priors. First, we have to choose the dimension of our latent space $X$. We choose $dim(X)=2$, since we would like our model to disentangle 'capture time' ($1$, $2$, $4$, $8$, $16$, $32$, and $64$) from cell branching types (TE, ICM, PE, EPI). Next, when we set the mean of our prior over $X$, we set the first dimension to be equal to the observed capture time. This will help the GPLVM discover the structure we are interested in and will make it more likely that that structure will be axis-aligned in a way that is easier for us to interpret."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "capture_time = y.new_tensor([int(cell_name.split(\" \")[0]) for cell_name in df.index.values])\n",
    "# we scale the time into the interval [0, 1]\n",
    "time = capture_time.log2() / 6\n",
    "\n",
    "# we setup the mean of our prior over X\n",
    "X_prior_mean = torch.zeros(y.size(1), 2)  # shape: 437 x 2\n",
    "X_prior_mean[:, 0] = time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use a sparse version of Gaussian process inference to make training faster. Remember that we also need to define $X$ as a `Parameter` so that we can set a prior and guide (variational distribution) for it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "kernel = gp.kernels.RBF(input_dim=2, lengthscale=torch.ones(2))\n",
    "\n",
    "# we clone here so that we don't change our prior during the course of training\n",
    "X = Parameter(X_prior_mean.clone())\n",
    "\n",
    "# we will use SparseGPRegression model with num_inducing=32;\n",
    "# initial values for Xu are sampled randomly from X_prior_mean\n",
    "Xu = stats.resample(X_prior_mean.clone(), 32)\n",
    "gplvm = gp.models.SparseGPRegression(X, y, kernel, Xu, noise=torch.tensor(0.01), jitter=1e-5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use the [autoguide()](http://docs.pyro.ai/en/dev/contrib.gp.html#pyro.contrib.gp.parameterized.Parameterized.autoguide) method from the [Parameterized](http://docs.pyro.ai/en/dev/contrib.gp.html#module-pyro.contrib.gp.parameterized) class to set an auto Normal guide for $X$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we use `.to_event()` to tell Pyro that the prior distribution for X has no batch_shape\n",
    "gplvm.X = pyro.nn.PyroSample(dist.Normal(X_prior_mean, 0.1).to_event())\n",
    "gplvm.autoguide(\"X\", dist.Normal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inference"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As mentioned in the [Gaussian Processes tutorial](gp.ipynb), we can use the helper function [gp.util.train](http://docs.pyro.ai/en/dev/contrib.gp.html#pyro.contrib.gp.util.train) to train a Pyro GP module. By default, this helper function uses the Adam optimizer with a learning rate of `0.01`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAD4CAYAAADCb7BPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAex0lEQVR4nO3dfZRcdZ3n8fe3qqs7z88NhCSYoFEIoBJiyKjDURghBMcwZ+AszAMZhzVzFGdUdmYMR1ccXWbRwcVhB3EiRAjj8iDjLCwGYhbw+EQCHXkM2ZgmiaRNSBryQEin093V3/3j/qpT3amurufbnfq8zqlTt373d+/v27dDf7gPda+5OyIiIsVKxF2AiIiMTAoQEREpiQJERERKogAREZGSKEBERKQkDXEXUCvTpk3z2bNnx12GiMiIsnHjxjfcvTnXvLoJkNmzZ9PS0hJ3GSIiI4qZ/XaweTqEJSIiJVGAiIhISRQgIiJSEgWIiIiURAEiIiIlUYCIiEhJFCAiIlISBcgQnt2xj2/9ZAvd6d64SxERGVYUIEN47rX9/M8nWznaowAREcmmABlCKhltoh7tgYiI9KMAGUImQLoUICIi/ShAhtAYAqQ7rUf/iohkU4AMoSFpAHTrHIiISD8KkCGk+vZAFCAiItkUIENI6RCWiEhOCpAhNDaEQ1jaAxER6UcBMgQdwhIRyW3IADGzVWa218xezmqbYmbrzGxreJ8c2s3MbjOzVjN70czmZy2zLPTfambLstrPM7OXwjK3mZmVOkY16DJeEZHcCtkDuRtYPKBtBfCEu88FngifAS4F5obXcuAOiMIAuBE4H1gI3JgJhNBnedZyi0sZo1pSmauwdA5ERKSfIQPE3X8G7BvQvBS4J0zfA1ye1b7aI+uBSWY2HbgEWOfu+9x9P7AOWBzmTXD3p93dgdUD1lXMGFXRdwhLl/GKiPRT6jmQk919N0B4Pym0zwB2ZvVrC2352ttytJcyxnHMbLmZtZhZS3t7e1E/YEbfrUx6FSAiItkqfRLdcrR5Ce2ljHF8o/tKd1/g7guam5uHWG1umQDRzRRFRPorNUD2ZA4bhfe9ob0NmJXVbyawa4j2mTnaSxmjKkalQoB0K0BERLKVGiCPAJkrqZYBD2e1XxOulFoEHAyHn9YCF5vZ5HDy/GJgbZh3yMwWhauvrhmwrmLGqIpRqSQAnT3pag0hIjIiNQzVwczuAz4CTDOzNqKrqW4GHjSza4HXgCtD9zXAEqAV6AA+CeDu+8zs68Czod/X3D1zYv7TRFd6jQYeCy+KHaNaRmcCpFsBIiKSbcgAcferB5l1UY6+Dlw3yHpWAatytLcAZ+dof7PYMaqhbw9Eh7BERPrRN9GHkExE5+wffbFqp1lEREYkBUiBfrPn7bhLEBEZVhQgBZoxaXTcJYiIDCtDngMROO8dk/su5xURkYj+KhZgVCqh74GIiAygAClAU0NS3wMRERlAAVIA7YGIiBxPAVIA7YGIiBxPAVKAiaNTHOjojrsMEZFhRQFSgDGNSY50aQ9ERCSbAqQATQ1JenqdHj3WVkSkjwKkAE0pPRddRGQgBUgBmhr0TBARkYEUIAVoaojuyKs9EBGRYxQgBdAeiIjI8RQgBcicAzmq74KIiPRRgBQg81TCDl3KKyLSRwFSgImjUwAcPKIvE4qIZChACjAhBMhbnQoQEZEMBUgBMoew9G10EZFjFCAFGBUCpLNHV2GJiGQoQAqQeRrh0W7tgYiIZChACjBKh7BERI6jAClAKpkgmTA9E0REJIsCpECjU0k69U10EZE+CpACjUolOKJzICIifRQgBWpqSNKpcyAiIn0UIAWaNCalb6KLiGRRgBRo6rgm3nj7aNxliIgMG2UFiJl9wcw2mdnLZnafmY0yszlmtsHMtprZA2bWGPo2hc+tYf7srPXcENq3mNklWe2LQ1urma3Ias85RjVNHdvIvo6uag8jIjJilBwgZjYD+BtggbufDSSBq4BvALe6+1xgP3BtWORaYL+7vwu4NfTDzOaF5c4CFgPfMbOkmSWB24FLgXnA1aEvecaomqaGBF36JrqISJ9yD2E1AKPNrAEYA+wGLgQeCvPvAS4P00vDZ8L8i8zMQvv97n7U3bcDrcDC8Gp1923u3gXcDywNyww2RtWkkgm6017tYURERoySA8TdfwfcArxGFBwHgY3AAXfvCd3agBlhegawMyzbE/pPzW4fsMxg7VPzjNGPmS03sxYza2lvby/1RwWgUXsgIiL9lHMIazLR3sMc4FRgLNHhpoEy/9tug8yrVPvxje4r3X2Buy9obm7O1aVgqWRCz0QXEclSziGsPwC2u3u7u3cDPwI+CEwKh7QAZgK7wnQbMAsgzJ8I7MtuH7DMYO1v5Bmjatr2d9DV00uPQkREBCgvQF4DFpnZmHBe4iLgFeAp4IrQZxnwcJh+JHwmzH/S3T20XxWu0poDzAWeAZ4F5oYrrhqJTrQ/EpYZbIyqefTF3QDsePNwtYcSERkRyjkHsoHoRPavgZfCulYCXwSuN7NWovMVd4VF7gKmhvbrgRVhPZuAB4nC53HgOndPh3McnwXWApuBB0Nf8oxRNd+84r1E9VZ7JBGRkaFh6C6Dc/cbgRsHNG8juoJqYN9O4MpB1nMTcFOO9jXAmhztOceopqljo6+aHNbtTEREAH0TvWBjGqOs7ejqGaKniEh9UIAUaExj9FCpjqPaAxERAQVIwcY2hQDRLd1FRAAFSMFGh0NYR3QIS0QEUIAUbEx4Lvret3RHXhERUIAUbEw4hPWtdb+JuRIRkeFBAVKgxqQ2lYhINv1VLFD0ZXsREcko64uE9eb05rHMnjo27jJERIYFBUgRtrUfZlu77oUlIgI6hFUUHcUSETlGAVKE8+dMibsEEZFhQwFShPXb9gHoyYQiIihAStLZo9uZiIgoQIrwnz88B4Cj3doDERFRgBRh9rToEl7XU6VERBQgxUgmosuw0goQEREFSDGS4TredK8CREREAVKEn7zyOgBrN+2JuRIRkfgpQIrwavgW+rPb98VciYhI/BQgRRgdngny+KbXY65ERCR+CpAiXHr2KXGXICIybChAivCpC06PuwQRkWFDd+MtwqhUkpmTR3PuaZPjLkVEJHbaAynSuKYGOrt1KxMREQVIkUalkgoQEREUIEXr6ull8+634i5DRCR2OgdSpFcUHiIigPZARESkRGUFiJlNMrOHzOz/mdlmM/s9M5tiZuvMbGt4nxz6mpndZmatZvaimc3PWs+y0H+rmS3Laj/PzF4Ky9xmFt2MarAxaqknrVu6i0h9K3cP5J+Bx939DOB9wGZgBfCEu88FngifAS4F5obXcuAOiMIAuBE4H1gI3JgVCHeEvpnlFof2wcaomd0HO2s9pIjIsFJygJjZBOAC4C4Ad+9y9wPAUuCe0O0e4PIwvRRY7ZH1wCQzmw5cAqxz933uvh9YBywO8ya4+9MePYBj9YB15Rqj6r7xx+cAcFSPtRWROlfOHsjpQDvwfTN7zszuNLOxwMnuvhsgvJ8U+s8AdmYt3xba8rW35Wgnzxj9mNlyM2sxs5b29vbSf9Is08Y1AXD4aE9F1iciMlKVEyANwHzgDnc/FzhM/kNJlqPNS2gvmLuvdPcF7r6gubm5mEUHNbYpunBNASIi9a6cAGkD2tx9Q/j8EFGg7AmHnwjve7P6z8pafiawa4j2mTnayTNG1Y0LAfK2AkRE6lzJAeLurwM7zew9oeki4BXgESBzJdUy4OEw/QhwTbgaaxFwMBx+WgtcbGaTw8nzi4G1Yd4hM1sUrr66ZsC6co1RdaNS0Sb7159tq9WQIiLDUrlfJPxr4Adm1ghsAz5JFEoPmtm1wGvAlaHvGmAJ0Ap0hL64+z4z+zrwbOj3NXfPPLHp08DdwGjgsfACuHmQMapu5uQxAGz87f5aDSkiMiyVFSDu/jywIMesi3L0deC6QdazCliVo70FODtH+5u5xqiFpgZ991JEBPRN9KKF7zKKiNQ93QurBO+dOZGpYxvjLkNEJFbaAylBYzJBl25lIiJ1TnsgJWjRCXQREe2BiIhIaRQgZejtLeqL8SIiJxQFSBme3vZm3CWIiMRGAVKCr/7hPCB6PrqISL1SgJTgXSeNByCtQ1giUscUICVIhK32xttH4y1ERCRGCpASvLDzIAB/98MXYq5ERCQ+CpASjBsVngnSlY65EhGR+ChASnDVB2YN3UlE5ASnAClBKhlttukTR8VciYhIfBQgZdh9sDPuEkREYqMAERGRkihAytTRpWeji0h9UoCU6TtPvRp3CSIisVCAlGnSmFTcJYiIxEIBUqL7PrUIgDnTxsZciYhIPBQgJRofvkzYo/thiUidUoCUqCFpgJ4JIiL1SwFSJj0TRETqlQKkRB52PFY//dt4CxERiYkCpERnTp8AwMXzTo65EhGReChAynD+nCkc6OiOuwwRkVg0xF3ASLZh+764SxARiY32QMrwnpOjR9t2p3tjrkREpPYUIGW4csFMQI+2FZH6VHaAmFnSzJ4zs0fD5zlmtsHMtprZA2bWGNqbwufWMH921jpuCO1bzOySrPbFoa3VzFZktecco9bWvLQbgG/95DdxDC8iEqtK7IF8Dtic9fkbwK3uPhfYD1wb2q8F9rv7u4BbQz/MbB5wFXAWsBj4TgilJHA7cCkwD7g69M03Rk0tOWc6ADMmjY5jeBGRWJUVIGY2E7gMuDN8NuBC4KHQ5R7g8jC9NHwmzL8o9F8K3O/uR919O9AKLAyvVnff5u5dwP3A0iHGqKnFZ58CwJuHdQhLROpPuXsg3wb+HsicRZ4KHHD3zEMy2oAZYXoGsBMgzD8Y+ve1D1hmsPZ8Y/RjZsvNrMXMWtrb20v9GQc1ZWx05Ozf1r9W8XWLiAx3JQeImX0c2OvuG7Obc3T1IeZVqv34RveV7r7A3Rc0Nzfn6lKWMY26ClpE6lc5fwE/BHzCzJYAo4AJRHskk8ysIewhzAR2hf5twCygzcwagInAvqz2jOxlcrW/kWeMmhvf1MCho3oqoYjUn5L3QNz9Bnef6e6ziU6CP+nufwo8BVwRui0DHg7Tj4TPhPlPuruH9qvCVVpzgLnAM8CzwNxwxVVjGOORsMxgY9Tc6MZkXEOLiMSqGt8D+SJwvZm1Ep2vuCu03wVMDe3XAysA3H0T8CDwCvA4cJ27p8PexWeBtURXeT0Y+uYbo+b+5PzTAEjrtu4iUmcqchDf3X8K/DRMbyO6gmpgn07gykGWvwm4KUf7GmBNjvacY8ThV63R7dxfaDvA/NMmx1yNiEjt6JvoZXpmR3Q/rL+6d+MQPUVETiwKkDJdvTA6hNV+SN8FEZH6ogAp0x+dm/MrKCIiJzwFSJnOPW1S3CWIiMRC34QrUyqZwAxGp3Q5r4jUF+2BVMCpE0fT0ZXmhZ0H4i5FRKRmFCAV8LsDRwC4auX6mCsREakdBUgFnD1jAnDsiiwRkXqgAKmA+z61CIBVv9wecyUiIrWjAKmAcU26FkFE6o8CpAKiZ1yJiNQXBYiIiJREAVJh2984HHcJIiI1oQCpsEOd3XGXICJSEwqQCmlsiDblz7e+EXMlIiK1oQCpkC9fdiYA/7R2S8yViIjUhgKkQs44ZULcJYiI1JQCpEIWzpnSN/3bN3UiXUROfAqQKvj+L3fEXYKISNUpQKrg7l/tiLsEEZGqU4BU0F9dcHrcJYiI1IwCpILOmD4+7hJERGpGAVJB80+bHHcJIiI1owCpoHdMHds3/cqut2KsRESk+hQgVbJljwJERE5sCpAKu/TsUwD4b49ujrkSEZHqUoBU2Jc/Pg+ANw93xVyJiEh1KUAqbNq4xr7pI13pGCsREakuBUiFNTUk+6b/9M71MVYiIlJdJQeImc0ys6fMbLOZbTKzz4X2KWa2zsy2hvfJod3M7DYzazWzF81sfta6loX+W81sWVb7eWb2UljmNgvPjh1sjOFi9V8uBODXrx2IuRIRkeopZw+kB/gv7n4msAi4zszmASuAJ9x9LvBE+AxwKTA3vJYDd0AUBsCNwPnAQuDGrEC4I/TNLLc4tA82xrBwwbub4y5BRKTqSg4Qd9/t7r8O04eAzcAMYClwT+h2D3B5mF4KrPbIemCSmU0HLgHWufs+d98PrAMWh3kT3P1pd3dg9YB15Rpj2Ln/mdfiLkFEpCoqcg7EzGYD5wIbgJPdfTdEIQOcFLrNAHZmLdYW2vK1t+VoJ88YA+tabmYtZtbS3t5e6o9XlhU/eimWcUVEqq3sADGzccC/A59393zfnrMcbV5Ce8HcfaW7L3D3Bc3NtT2s9NJXL67peCIitVZWgJhZiig8fuDuPwrNe8LhJ8L73tDeBszKWnwmsGuI9pk52vONMWyMH5Xqm/7Vq3pOuoiceMq5CsuAu4DN7v4/smY9AmSupFoGPJzVfk24GmsRcDAcfloLXGxmk8PJ84uBtWHeITNbFMa6ZsC6co0xrPzX8KXCP/nehpgrERGpvHL2QD4E/DlwoZk9H15LgJuBj5nZVuBj4TPAGmAb0Ap8D/gMgLvvA74OPBteXwttAJ8G7gzLvAo8FtoHG2NY+eQHZ/dNay9ERE40DaUu6O6/IPd5CoCLcvR34LpB1rUKWJWjvQU4O0f7m7nGGG4SiWOb5/5ndvLBd06LsRoRkcrSN9Gr7DMfeScAj7ywa4ieIiIjiwKkyj7/B+/umz7rK4/HWImISGUpQKqssSHB314chcjhrjQ96d6YKxIRqQwFSA1cfNYpfdPXP/hCjJWIiFSOAqQG3n3yeM44ZTwQnQvZffBIzBWJiJRPAVIjj3/+gr7p3/vvT8ZYiYhIZShAaujCM47dsmv2ih/HWImISPkUIDW06i8+wBXnHbs7y5Xf/VWM1YiIlEcBUmO3XPm+vulnd+xnw7Y3Y6xGRKR0CpAY/J/Pfrhv+j+tXM9VK5+OsRoRkdIoQGJwzsyJvPqPS/o+r9+2j1vWbomxIhGR4ilAYpJMGDtuvqzv87881crS239JdMswEZHhTwESs603XcqZ0ycA8MLOA8y5YQ3d+ra6iIwACpCYpZIJHvvc7/e7Omvulx7TzRdFZNhTgAwTt1z5Pn5/7rHbvf/Nfc8xe8WP6ejqibEqEZHBKUCGkXuvPZ9ffPGj/drmfWUtX/7fL8VUkYjI4KxeTtouWLDAW1pa4i6jYP/xXBtfeKD/jRfHNiZ5+R8uIXrCr4hI9ZnZRndfkGue9kCGqT86dybb/nEJ5542qa/tcFeaOTesYfaKH7PqF9t1sl1EYqU9kBGgszvNLWu3cOcvtuec/90/m88lZ52iPRMRqbh8eyAKkBHm4JFuPnDT/6WrJ/fex5JzTmHF4jM5beqYGlcmIiciBQgnToBk+8mm11l+78a8fa5eOIu/u+QMpoxtrFFVInIiUYBwYgZItlfb3+ZARzd/fMfQd/j9zEfeyUfPOIkPzJ5Sg8pEZCRTgHDiB8hAO/d18Nf3PcfzOw8UvMxffHA2V5w3k3edNI5RqWQVqxORkUIBQv0FSC4HOrr41OoWnt2xv6Tll77/VP5s0Ts4Z8ZEGhJGQ1IX8Ymc6BQgKEDy2f7GYX6+tZ2te97m3vW/LXr5j76nmVlTxtDT6+w52Mkn3n8qTQ1J1m56nQvPOIlzT5vEtHFNNCYTmKGrxURGEAUICpBSHT7aw64DR3js5dd5Zvs+enp7OXikh8273+rrM21cEx1dPXR0pYte/0ff08zE0SlGNyZ5+Pld/Pmid7DrYCfjmpLMO3UizeOaSCWNpoYke97q5JyZEznSleakCU2MH5WiccBeUGOD9opEKilfgDTUuhgZWcY2NTD35PHMPXn8kH27070c6Ohm76FOdrzRwXd+2soru99i8Vmn8PS2NznQ0X3cMrsOdLLl9UO8fTQKoH/92basuTtLqnnauEYaEgmSCSOVjA61RYfcjIbEselUMuqTqy0V+h63jkQirCdqSyUt6hP6NiQNM8PdSSUTpJIJEmGHK2GGWfQevaK9sYRBIpH9ObTl6581LyOzd5c0I5EAI+pjRO3pXscMUokoZC0B7mSNcWx9mXfHMezY+oG0O6lEgp6wvmi84/cq3V17myc4BYhUTCqZoHl8E83jmzjr1Ilc9t7pRS3f2+sc6U7zVmc3PWnHHdrf7iTdC73uPPfaASaNSfHca/t536xJHDzSTfYO9M59HTQkjV6HdNrp7u2lJ+2ke53udC89vR690lH74Z6e8NnpCX0z87t7s5bLrKO3lzrZYS9awiB70yTNSLuTNOtrz0RJJlP6BZP1/5zpnwkgy1pBdvtgBs62fvNskPbj1pJzXr5lrIBlBtbQv728dQ9cb+bT5z/2bj7xvlNzjlkOBYgMG4mEMbapgbFNx/5ZZn8hctHpUwG4euFpNa8tI92bFTaZ4MkKpu60A86hzh6aGpKke520OwmDXo+C0N2j6d7ove+ze5ifmR7Q37P7O70hWPv+cHu0x9AT1kvWvN5eJ5lM4O50px0jWrZv0RxjwrE/SB7qcqKHoXWne0ma9fXP7N1YVAbpXidhRq9H75lKMkN61rieVWdm3MxYx6b7L5/pm+sP8cDD8t5vXna752zPtwz5lilz3dnLkHcZH6R98GUmj0lRDSM2QMxsMfDPQBK4091vjrkkqQPJhJFMJGkasf/liFTOiDzjaGZJ4HbgUmAecLWZzYu3KhGR+jIiAwRYCLS6+zZ37wLuB5bGXJOISF0ZqQEyg/6X6LSFtn7MbLmZtZhZS3t7e82KExGpByM1QHJdwnDc9THuvtLdF7j7gubm5hqUJSJSP0ZqgLQBs7I+zwR2xVSLiEhdGqkB8iww18zmmFkjcBXwSMw1iYjUlRF5MaK795jZZ4G1RJfxrnL3TTGXJSJSV0ZkgAC4+xpgTdx1iIjUq7q5maKZtQPF32o2Mg14o4LlVIrqKs5wrQuGb22qqzgnYl3vcPecVyHVTYCUw8xaBrsbZZxUV3GGa10wfGtTXcWpt7pG6kl0ERGJmQJERERKogApzMq4CxiE6irOcK0Lhm9tqqs4dVWXzoGIiEhJtAciIiIlUYCIiEhJFCBDMLPFZrbFzFrNbEUM4+8ws5fM7HkzawltU8xsnZltDe+TQ7uZ2W2h1hfNbH4F61hlZnvN7OWstqLrMLNlof9WM1tWpbq+ama/C9vseTNbkjXvhlDXFjO7JKu9or9nM5tlZk+Z2WYz22RmnwvtsW6zPHXFus3MbJSZPWNmL4S6/iG0zzGzDeFnfyDcuggzawqfW8P82UPVW+G67jaz7Vnb6/2hvWb/9sM6k2b2nJk9Gj7Xdnt5eHylXse/iG6T8ipwOtAIvADMq3ENO4BpA9q+CawI0yuAb4TpJcBjRHcrXgRsqGAdFwDzgZdLrQOYAmwL75PD9OQq1PVV4G9z9J0XfodNwJzwu01W4/cMTAfmh+nxwG/C+LFuszx1xbrNws89LkyngA1hOzwIXBXavwt8Okx/BvhumL4KeCBfvVWo627gihz9a/ZvP6z3euB/AY+GzzXdXtoDyW+4PrhqKXBPmL4HuDyrfbVH1gOTzGx6JQZ0958B+8qs4xJgnbvvc/f9wDpgcRXqGsxS4H53P+ru24FWot9xxX/P7r7b3X8dpg8Bm4meWRPrNstT12Bqss3Cz/12+JgKLwcuBB4K7QO3V2Y7PgRcZGaWp95K1zWYmv3bN7OZwGXAneGzUePtpQDJr6AHV1WZAz8xs41mtjy0nezuuyH6gwCcFNprXW+xddSyvs+GQwirMoeJ4qorHC44l+j/XofNNhtQF8S8zcLhmOeBvUR/YF8FDrh7T44x+sYP8w8CU2tRl7tnttdNYXvdamZNA+saMH41fo/fBv4e6A2fp1Lj7aUAya+gB1dV2YfcfT7R89+vM7ML8vQdDvXC4HXUqr47gHcC7wd2A9+Kqy4zGwf8O/B5d38rX9da1pajrti3mbun3f39RM/3WQicmWeM2Ooys7OBG4AzgA8QHZb6Yi3rMrOPA3vdfWN2c54xqlKXAiS/2B9c5e67wvte4D+I/sPakzk0Fd73hu61rrfYOmpSn7vvCf/R9wLf49gueU3rMrMU0R/pH7j7j0Jz7NssV13DZZuFWg4APyU6hzDJzDJ3Dc8eo2/8MH8i0aHMWtS1OBwKdHc/Cnyf2m+vDwGfMLMdRIcPLyTaI6nt9ir3JM6J/CK63f02opNLmROFZ9Vw/LHA+KzpXxEdN/0n+p+I/WaYvoz+J/CeqXA9s+l/srqoOoj+T2070UnEyWF6ShXqmp41/QWiY7wAZ9H/hOE2opPBFf89h599NfDtAe2xbrM8dcW6zYBmYFKYHg38HPg48EP6nxT+TJi+jv4nhR/MV28V6pqetT2/Ddwcx7/9sO6PcOwkek23V8X+uJyoL6KrKn5DdDz2SzUe+/Twy30B2JQZn+jY5RPA1vA+JbQbcHuo9SVgQQVruY/o0EY30f+1XFtKHcBfEp2oawU+WaW67g3jvkj0pMrsP45fCnVtAS6t1u8Z+DDRoYAXgefDa0nc2yxPXbFuM+C9wHNh/JeBr2T9N/BM+Nl/CDSF9lHhc2uYf/pQ9Va4rifD9noZ+DeOXalVs3/7Wev9CMcCpKbbS7cyERGRkugciIiIlEQBIiIiJVGAiIhISRQgIiJSEgWIiIiURAEiIiIlUYCIiEhJ/j94SQJDsbLgXwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# note that training is expected to take a minute or so\n",
    "losses = gp.util.train(gplvm, num_steps=4000)\n",
    "\n",
    "# let's plot the loss curve after 4000 steps of training\n",
    "plt.plot(losses)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After inference, the mean and standard deviation of the approximated posterior $q(X) \\sim p(X | y)$ will be stored in the parameters `X_loc` and `X_scale`. To get a sample from $q(X)$, we need to set the `mode` of `gplvm` to `\"guide\"`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "gplvm.mode = \"guide\"\n",
    "X = gplvm.X  # draw a sample from the guide of the variable X"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizing the result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let’s see what we got by applying GPLVM to our dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf4AAAGJCAYAAABrSFFcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOydeZgU1dX/P3f2YYBBEGQUlMUFFEcQBE0cREdBowRDfA1KXhdieI36hGj0527QxIhGg6hk4TUajInE1w1xAeIAEYkoKDqibIKoyCD7sMzCLPf3R3XNdFdXVVf13jPn8zz9zPStqlu3l5lz77nfc47SWiMIgiAIQvsgK9UDEARBEAQheYjhFwRBEIR2hBh+QRAEQWhHiOEXBEEQhHaEGH5BEARBaEeI4RcEQRCEdoQYfiHhKKXOUErNUUptUUodUkrtU0qtUEr9WilVYjlXBz0alVJfKKWeVkr1CjrnqsDxYx3ut1opVekynhMC198VeP6bwPODSqmONudfEzSmPtG+D/FEKTVeKbVUKbVDKVWjlNqslHpZKTU66Bxz3L3c+opxHAm5h1LqGKXUTKXU50qpOqXUAaXU+0qp25VSnX32dW5gjGcGtb2jlHorjuN9x/LdrVJKvaGUOs3m3IivTSmVY+mvOfBZv6yUGhjDOI8N9PfjKK69SSl1cbT3FtIHMfxCQlFK/RJYBnQH7gLOBSYAC4DJwFM2l/0VOAMYBTwCfB+oUEoVerztM8DJSqlTHI5fAWjgb5b2JuCHDufv93jvhKOUugl4EVgLXA2MBe7H+HseFXTqXIz3cXuShxgTSqlRQCVQDkwHzsf4XF4Hfg7cnbLBubMK4/3+DnATcAzwtlLqBPOEKF7bXwJ9jgSmAmXAfL+TnzhxEyCGvy2gtZaHPBLyAM4GmoHpDseLgKssbRr4jaXtykD7+MDzqwLPj3Xo90igEXjE5pgCNgOLgtp+E+jvr8BblvP7BF7D04Fz+qTB+7oV+D+HY1lJHss1gfelV5z6OxzYCbwDdLA53gk412ef5wbGeGZQ2zvWzzrGcb8DLLG09Q/c9/d+XxuQE7h2quUc82/hkijHeWzg+h9Hce0W4K/J/H7JIzEPWfELieRWjH90t9od1Fof1Fr/1UM/KwI/bV37Nv1uBd4CLldKZVsOj8JYiT1jc+kzwNlKqaOC2q4ANmF4LTyhlLpSKVWplKoPuGdnK6WOsJyzRSn1V6XURKXU2sA2wwql1Hc83KIrsM3ugNa6OegeYW54P/cNuHa/VErVKqWWK6VGBK5/0sN7cG3gPagLvAf/q5Tq4uG1TQa6ATdorWtsXt9+rXWLi14p1VEp9bvAVschpdQmpdRtSinl4V4RUUr1V0q9GdhO2a6U+r1S6jov2xta643AHlq/t75emwMfBn4e7WHsRUqpPymldge2E17BmBRbzxuhlHox8NnWKqXWKWP7qyDonC3AUcCVQdsPTwaOHa+UejbwGdQqpTYGtjK8fN5CCshJ9QCEtolSKgc4C3hJa30oxu76Bn7u9XHNbOAfGKu9BUHt/w3UAC/YXLMYY1UzEXgo6HzrloAjSqnrgJmBe98K9AZ+CwxXSg21/MM/GxgI3AkcwvA8vKaU6qO13udym/eBSUqpzcCrWusNXsfn9b5KqWsxtln+F2Nb4Vjgn0BEF7NS6mFgCvAocDPQC2Mr4iSl1JnBkxMbzgW+1lp/5OE+ucBC4Hjg18BqDDf7vcBhOEw4vaKUygcqgFzgZxiT2OuA//J4/WFAF1q/t55fmwt9Aj83ejj3SWA8xhbBB8AY4Fmb847BmFA8DRwATgLuCdzL1AKMxfg7WoHxXkPrFtJRwJfA/9E60bkDOAVo0VUIaUSqXQ7yaJsP4AgMl+IDNsdygh+WYxrDSOQABcDpwBrgIHBk4JyrcHH1B84pxPiH+3dL2z7gGcu5vzH+FDQYRvqTwO/fwXDz96PVpd3H5Z45wA7CtwtGBa69LqhtC7ALKA5qOz1w3qUR3tsBGEZOBx47CExyLOeFueG93BfIxthOeNXS36WB8550ugeGe7sJuMNy7VmB8y6K8No2AEs9fseuDvT5HUv7r4B6oFvgeVSufgxjr4HTgtqyMbQV1vf1HeDfge9ALobxmxf8mn2+NtPVfx+tfwvDgc8C98qJcP2Jge/uzZb2/8XF1Y+xFZaD8TfWBHSxfHf+6nHs5nf+ZC+vVx7JfYirX0gUtq5WpVRPoCH4EfAOBHNH4Fgt8G7g9+9pw4XvCa11LcYK5GKlVKdA8w8w9lHt3PwmzwCDlFJDMNz8y7TWmzze9kSMfdyQVZXWegnwDYbxC2aZ1ro66PkngZ+ublyt9VqM1dQojIlKJYZA7F9Kqds8jDPSfY8BSjDev2BexjAmbozGEBn+PaBMzwl8vsswPC0jAZRS2cHHo3TNn4+x8n3fcq+FQB4wIoo+gzkD+EJrbW41obVuIvx9MRmJ8V09hGHkhwM/1Vq/FsMY7qb1b+E9jAnAxVrrxgjXnY7xN/i8pX2O9USlVJfAdskmjAlTA8bqPwsP22tKqXyl1F2BraPawPWLA4dPcLlUSBFi+IVEsROoI9yI7QROCzz+1+HapwLHhwCHa61Ltdb/jmIMs4EOtCr1r8BYtSxyuiBgVFdgrGR/hPskwUrXwM8qm2Pbgo6b7LY8rw/8LCACWusmrfW/tdZ3aq3LMbwSnwL3eVB8R7qvGWIZEg2gtW6wudZKj8DPzVgmeBifRbfA8X9bjt0ZaP8aY+LhhR4YHgbrff4TON7N4TqvlADf2rTbtYHhLj8NGIaxvdFTax2sh/Dz2kz+N9BnGcbq/xjg7x6uMz9D61jtxj4b+CnG1sx5gfv9PHAs4ncRY1vsHoy/lQsxJjzmdoiX64UkI3v8QkLQWjcqpd4GzlNK5enAPn9gpbISQCl1kcPlVVrrlXEYwztKqY3Afyul5mO4fB/S7nvMYPwDm4GxcrOumNwwjWJPm2M9gc999OULrfU3SqmnMPblj6VVBBYN5sSlR3BjYE/dOnmxsivwsxxjW8XKzsDPn2B4X0y+Cfx8C0NgOVhH3gvfhfGeXuZw/IsI10eiCvs96iNs2gD2R/je+nltJluD+nwnIFa9Uyn1A631yy7XmZ/hEcBXQe1WkWkH4CLgTq31Y0HtQzyOD4zw3Ke01r8Nul6EfWmMrPiFRPIQhuv7wRSO4RkMl/itGPuzXlbwzwGvYugTqiOdHMRnGIZtQnCjUuosDAFUNF6LMJRSvR0ODQj8tFX8++BLDMNhFbGNJ/L/jIUYe7u9tdYrbR6bAbTW6yztpqGahWHQn1A2eRsCSvVzA0/nY3iUqh3utct6vU/eBfoqpYYF3T8bj+I+G/y8Nid+i7Fq/1WE7ZHlBHQblvYJlucFGJ9pQ9A4FMYev5V6DJ2MlcLg6wNc7TI2IcXIil9IGFrrisCe8zSlVCmG0f0C45/N8Rj/hA5i/IOKhvOVUlYjV621/lfQ82cwVM1TgPcDrvxI496FoQfwRcDL8StgplJqNsYEohfGP+u1GC7VeLBGKbUQeAXDpd4ZY9X2U+AffrQQdmitm5RS9wF/VErNwoiAOA74fxiJjBw9Jlrr9QFV/x+VkWHubQyD0Rtj//+PWuulLtfvVEpdgpF86EOl1BMYQsZ8DBfyzzCEjG9hfLZXAYsD9/wEY2//WIykTxdprevDbuKdpzAmjK8ope6kVdVfFE1nPl+bUx81SqkHMNzy3w/0ZXfeZ0qpfwL3B3QPH2BoIsZYztutlFoJ/D+l1LcYqvxrsPdqfAacpZS6EGPysUNr/SWG2n+SUuozDM3FfwVej5CupFpdKI+2/wC+i+Ey/wbDfb4PYx/9XqDEcm5YAh+b/q6iVdFufay2OX9J4Nj1Dv21qPpd7hlR1R907pUYgrt6DGMxGzjCck6YQppWJfddEfq/HkMx/iWGjuIghmv/FiDXZsxWVb+n+wK/xHAT1wU+r+8GPrvfud0j6D14D0PQtx/DaDxOIDLDw3vYB/gDhiGpxwgzex9j8tEp6LxCjL3vdYHzdgXu+ysCyYyIIYEPxiRiPoa4bjvwewzjb6fqXxKv14ZDAp/AsfzA57Iywn2KgD9jGPMDGJOEkVhU/Rj6kPmBz2k78BjGpML6np0YeJ01BEV3YGTlfB4jimYPRvjrCOt95JE+DxX44ARBEFxRSp2BIZy7XGv9XKrHkyqUUtdgiO56a623pHo8guAXcfULghCGUqo/cC2wFGMleBJwO4aYzk1UJghCmiOGXxAEO2qBUoxtlS4YLtyFwG1a67oUjksQhBgRV78gCIIgtCMknE8QBEEQ2hFi+AVBEAShHdHm9/gPP/xw3adPn1QPQxAEQRCSxgcffLBTa93d7libN/x9+vRh5cqYs78KgiAIQsaglPrS6Zi4+gVBEAShHSGGXxAEQRDaEWL4BUEQBKEd0eb3+AVBEIS2QUNDA1u2bKGuTnJImRQUFNCrVy9yc3M9XyOGXxAEQcgItmzZQqdOnejTpw/uVYnbB1prdu3axZYtW+jbt6/n68TVLwiCIGQEdXV1dOvWTYx+AKUU3bp18+0BEcMvCIIgZAxi9EOJ5v0Qwy8IgiAIHpk0aRI9evRg0KBBqR5K1IjhFwRBEASPXHXVVcyfPz/Vw4gJMfyCIAhCm+TFbbsZ9p9PKVn8EcP+8ykvbtsdc58jR46ka9eucRhd6hBVvyAIgtDmeHHbbm5e9zW1zUbp+S31Ddy87msAftgzsw13rMiKXxAEwUrl8zB9EEztYvysfD7VIxJ88sCmqhajb1LbrHlgU1WKRpQ+yIpfEAQhmMrnYd7PoaHWeF79tfEcoPTS1I1L8MU39Q2+2tsTsuIXBEEIpuK+VqNv0lBrtAsZw1H59pnsnNrbE2L4BUEQgqne4q9dSEtu71dCYVZojHthluL2fiUx9XvZZZdxxhlnsG7dOnr16sVf/vKXmPpLBeLqFwRBCKa4l+Het2sXMgZTwPfApiq+qW/gqPxcbu9XErOw77nnnovH8FKKGH5BEIRgyu8J3eMHyC002oWM4oc9u7Z7Bb8d4uoXBEEIpvRSGPsYFPcGlPFz7GMi7BPaDLLiFwRBsFJ6qRh6oc0iK35BEARBaEeI4RcEQRCEdoQYfkEQBEFoR4jhFwRBEAQPfP3115x99tkMHDiQk046iRkzZqR6SFEh4j5BEARB8EBOTg6PPPIIp556Kvv372fo0KGcd955nHjiiakemi/SasWvlDpfKbVOKfW5Uuo2m+PXKqU+UUp9pJR6RymVWe+2IAiCkDSqts1l2bIyKhYdy7JlZVRtmxtTfyUlJZx66qkAdOrUiYEDB/LNN9/EY6hJJW0Mv1IqG5gJXACcCFxmY9j/obU+WWs9GHgI+H2ShykIgiBkAFXb5rJ27Z3U1W8FNHX1W1m79s6Yjb/J5s2bWbVqFSNGjIhLf8kkbQw/MBz4XGu9SWt9CJgDjAs+QWu9L+hpERBac1EQBEEQgE0bH6a5ObTYUnNzLZs2Phxz3wcOHOCHP/whjz76KJ07d465v2STTnv8RwHBCbK3AGFTKaXU9cBNQB5wjl1HSqnJwGSAo48+Ou4DFQRBENKbuvoqX+1eaWho4Ic//CETJ05k/PjxMfWVKtJpxa9s2sJW9FrrmVrr/sCtwF12HWmtZ2mth2mth3Xv3j3OwxQEQRDSnYJ8+yp8Tu1e0Frzk5/8hIEDB3LTTTdF3U+qSSfDvwXoHfS8F7DV5fw5wMUJHZEgCIKQkfTrfzNZWYUhbVlZhfTrf3PUfS5btoy//e1vLFq0iMGDBzN48GDeeOONWIeadNLJ1b8COE4p1Rf4BpgAXB58glLqOK31hsDTC4ENCIIgCIKFkp6GRGzTxoepq6+iIL+Efv1vbmmPhjPPPBOtM19aljaGX2vdqJS6AVgAZANPaa0/VUrdB6zUWr8K3KCUOhdoAPYAV6ZuxIIgCEI6U9JzXEyGvq2SNoYfQGv9BvCGpe2eoN+nJH1QgiAIgtCGSKc9fkEQBEEQEowYfkEQBEFoR4jhFwRBEIR2hBh+QRCETKHyeZg+CKZ2MX5WPp/qEQkZiBh+QRCETKDyeZj3c6j+GtDGz3k/F+OfApqamhgyZAgXXXRRqocSFWml6hcEQRACVD4PFfdB9RYo7gWHDkJDaO55GmqNc0ovTc0Y2ykzZsxg4MCB7Nu3L/LJaYgYfkEQhHTDXN2bhr76a+dzq7ckZ0wZSGVlJRUVFVRXV1NcXEx5eTmlpaUx9bllyxZef/117rzzTn7/+8wsECuufkEQhHSj4r7w1b0Txb0SO5YMpbKyknnz5lFdXQ1AdXU18+bNo7KyMqZ+f/GLX/DQQw+RlZW55jNzRy4IgtBW8bqKzy2E8nsin9cOqaiooKGhIaStoaGBioqKqPt87bXX6NGjB0OHDo11eClFDL8gCEK64bSKL+wKxb0BZfwc+5js7ztgrvS9tnth2bJlvPrqq/Tp04cJEyawaNEifvzjH0fdX6oQwy8IgpBulN9jrOaDyS2ECx6EG1fD1L3GTzH6jhQXF/tq98IDDzzAli1b2Lx5M3PmzOGcc87h2Wefjbq/VCGGXxAEId0ovdRYzcvqPmrKy8vJzc0NacvNzaW8vDxFI0ofRNUvCIKQaqyhe+X3GEY+1YbeaVwZgKnej7eq32TUqFGMGjUqLn0lGzH8giAIqcQudG/ez43fU2lk03VcPigtLY2boW9LiKtfEAQhldiF7pmJeZKFXSrgdBiXkBBkxS8I7Yg1SxezdM4z7N+1k07dDqdswhUMLDs71cNq3ziF7iUrMY/Tyt4pj4AkDMp4ZMUvCO2ENUsXs3DWE+zfuQO0Zv/OHSyc9QRrli5O9dDaN06he8lKzOO0slfZ9udLwqCMRwy/ILQTls55hsZD9SFtjYfqWTrnmRSNSACcQ/eSlZjHaQWvm1I7LiFhiOEXhHbC/l07fbULSSLW0L1YS/U6ehx6S0hhG0X2+AWhndCp2+GGm9+m3QnRBCSJaEP34qG8L78nfE/fXNmnQ0hhmjF9+nSefPJJlFKcfPLJPP300xQUFKR6WL6QFb8gtBPKJlxBTl5+SFtOXj5lE66wPV80ARlAPJT3mZAsyPRq7P0Kvv0UananZBjffPMNjz32GCtXrmT16tU0NTUxZ86clIwlFmTFLwjtBHOl7nUF76YJkFV/mhCviIBIK/tUJvKxejWaDrWWKe7Q1fXSg6u2s2/BZpr21pPdJZ/OY/pQNKRHTMNpbGyktraW3NxcampqOPLII2PqLxWI4ReEdsTAsrM9G23RBGQAxb1ajaC1PV6kOpGPnVdDN8P+KlfDf3DVdva+tAHd0AxA09569r60ASBq43/UUUdx8803c/TRR1NYWMjo0aMZPXp0VH2lEnH1C4Jgi9Pev5smQEgyyYgISHUiHyfvRdMh18v2LdjcYvRNdEMz+xZsjnooe/bsYe7cuXzxxRds3bqVgwcPSpEeQRDaDn41AUIKSMb+fKoTDDl5L7LzXC9r2lvvq90Lb731Fn379qV79+7k5uYyfvx4/vOf/0TdX6oQV78gCLb41QQIKSLRyvtkbCe4YRd1oLKgU4nrZdld8m2NfHaXfJuzvXH00UezfPlyampqKCwspKKigmHDhkXdX6oQwy8IgiN+NAFCG8Ut3C8ZmJMac2shO88w+hGEfZ3H9AnZ4wdQuVl0HtMn6qGMGDGCSy65hFNPPZWcnByGDBnC5MmTo+4vVSitdarH0IJS6nxgBpANPKm1nmY5fhNwDdAI7AAmaa2/dOtz2LBheuXKlQkasSAIggPRKOHTtQxumoxrzZo1DBw40PP5iVD1pyN274tS6gOtta07Im1W/EqpbGAmcB6wBVihlHpVa/1Z0GmrgGFa6xql1M+Ah4AfJX+0giAILkSjhE+1et6NDE3kUzSkR5s09LGSNoYfGA58rrXeBKCUmgOMA1oMv9Y6OHPIcuDHSR2hEBcqKyupqKigurqa4uJiysvLpWa20LZwU8I7GdBorhGEKEgnw38UEKwg2QKMcDn/J8CbCR2REHcqKyuZN28eDQ0NAFRXVzNv3jwAMf4JQtLupoBolPCpVs8L7YZ0CudTNm22AgSl1I+BYcDvHI5PVkqtVEqt3LEjPDe5kDoqKipajL5JQ0MDFRUVKRpR20bS7qaIaErtpro8r1diLQokpJx0MvxbgN5Bz3sBW60nKaXOBe4Evq+1tg3I1FrP0loP01oP6969e0IGK0RHdXW1r3YhNqQUb4qIJrFOqsvzOhFs6B/sC69cFwjv0606BDH+GUU6Gf4VwHFKqb5KqTxgAvBq8AlKqSHAnzGM/vYUjFGIkeLiYl/tQmw4pt3duYNHJoxl1vVXy+o/EUSTWMfLNclebZuCQ9PQ1+6G5lCPXVKz+AlxIW32+LXWjUqpG4AFGOF8T2mtP1VK3Qes1Fq/iuHa7wj8n1IK4Cut9fdTNmjBN+Xl5SF7/AC5ubmUl5encFRtF6dSvECI6x+Qff9440UJbxcmd+Nq53OTrfq3Exza0Y50CJMmTeK1116jR48erF7d+lk9/vjjPPHEE+Tk5HDhhRfy0EMPpXCU7qTTih+t9Rta6+O11v211vcH2u4JGH201udqrY/QWg8OPMToZxilpaWMHTu2ZYVfXFzM2LFjRdiXIOzS7loR13+KsK6mI7nNU5Ez36tBTzcdQgK56qqrmD9/fkjb4sWLmTt3LpWVlXz66afcfPPNKRqdN9JmxS+0H0pLS8XQJwlr2l0cEnYFbwn4jQKQqIEo8Ru+lwrVv1O63mDSQYfgQCK+myNHjmTz5s0hbX/84x+57bbbyM83Jtk9eqR37oC0WvELghB/BpadzeSZT/PLOfPodLi92NWsuOc3CkCiBjxitzfv15CnQvVvJzjMzoPCriSsKFCcSOZ3c/369SxdupQRI0Zw1llnsWLFirjfI56I4ReEdoSd619lZ9NQX8cjE8by5h+m+4oCkKgBDzi59AsPsz/fyZBHUv0nQvhnJzgcNxNu/QKm7jX0CGlo9CG5383Gxkb27NnD8uXL+d3vfsell15KOqXDtyKufkFoR1hd//lFHWmoq6Vu/34Ax39WbtEBftrbJU4u/ZxCw3B7LX4TXKzGmjPfr/DPT+79DE3X6/iddWiPhV69ejF+/HiUUgwfPpysrCx27txJuoaTi+EXMhpJ/+uf4Ip7s66/mvoD+yNeY24FWFFZWejmZtt2IYCT6752D4yf5a/4jZMR9qMXcJokfLUcPn3ZCNkDw51/wYPutQXSoHCPE04RLU7f5Vi4+OKLWbRoEaNGjWL9+vUcOnSIww93uU/NbthfBU2HPFcbjCdi+IWMRdL/thKtiMnr6qffkNNs72Fn9AHH9naJWz37SKtpr8bVj17AaZKw8i+hbbW7jWQ94H3yYHduiiibcAULZz0R4u7PycunbMIVMfV72WWXsWTJEnbu3EmvXr249957mTRpEpMmTWLQoEHk5eUxe/ZsAiHn4dTsNt4vHfgbaTrU+v1IkvEXwy9kLG7pf9uT4TdFTOY/uEix+cEGXCnlaS+yctECPlm8kObGxpB75HfsZOsxcBIRtkuirWfvx7i6TS6s+IkCaG6w9xpkQEEh67ZWvFT9zz33nG37s88+662D/VWtRt9ENxvtYviFtkq83POS/tfATcRk/Se3Zuli5v9pRosB9ypA0k1NYYUzGg/Vk5OfR05eftxXVW0Kt715N/wYVz+TCy8hesFYJwqVzztfn2aJfIK3tdKGpkP+2hOAbMQJScV0z5vG2XTPV1ZW+u5L0v8a+BExLZo9q8Xoh6AUKOV7pV534ACjJ99gXBe4fvTkG9Lvn22qKb3UUMD7UcL7DffLCVL8F3Z1DrOziw6wrZEWINhrYHohvJwbjBT2aSU7z197ApAVv5BU4umel/S/Bn5ETKZ6Pwyt+eU/XwPg95d93/Mefaduh7esqswthDdm/p6Kv85CKWNiIEl9gvAjiHNamassox/zOuuWAECjS5pdOw/EcaNh1d/CV51ZuaFeA7cUvnYehsrn4c1bWwWDkJZ6gKTSqSR0jx+Mz7RTSdKGICt+IanE0z0v6X8N7GLzY3G3ezX65j3WLF3MzGsu440nHmlJllJ/YL8xyWivSX3sVrh+U/TarswB3RR6XTxS+R59uhGfXxi0x1zYFS7+Q6hxdnPlm/c0x2W+3mCjH+342hIduhr5EMwVfnae8VxU/UJbpbi42NbIR+uel/S/3kRM5mrcifyOnVp+73R4d1sPQn7HTuQVFITcAwhTTtvhpDlokziJ8nIK/QnizLaXrzWMvdN1XrcEWrwNX2O49nXo+MY+ZiTmcSOSPiB4NR+pwI85vjQPC0wIHbom1dBbEcMvJBVxzycGNxGTVfVvRWVnU37V5JbnTmFQ5VdNDrvHrOuvjmj0TRKROCUtcVqBOxnB6q/htZtgw8Jw41d6Kbw02eG6gOH0ougP2w6wSDW9KvLtRIRWGmoD7v097n0V98qIsMC2iBh+IamYq/N0T7rTlhID2an+TTod3j3MO+AnDMqPMU9E4pS0JBple3AMvdX4RTLsXhT9Xsrrehm3VR8QFusRoHa3sVVg5+YPHl8GhAUGU1dXx8iRI6mvr6exsZFLLrmEe++9F4CJEyeycuVKcnNzGT58OH/+85/Jzc1tuXbBggXceuutAHz++eccddRRFBYWUlpayqRJkxg3bhx9+/ZtOf/hhx/m3HPPTcjrEMMvJJ10d8+3tcRAjsZZKSbPfNr2kNcwKCdhoZV2FeLnZKgLuxqiOy/17YONXyTD7iVc0ItR91rsJzjp0PRB7q7/7Dz7MLVTLvfmzUgz8vPzWbRoER07dqShoYEzzzyTCy64gNNPP52JEye2xPJffvnlPPnkk/zsZz9ruXbMmDGMGTMGgFGjRvHwww8zbNgwAJYsWUJZWRmvvfZaUl6HGH5BsNAWEgN5SdITjxW43bYAQHZ+Prl5ee1C1b/+vW28O3cjB3bX07FrPqOHTaFkzd3hhvqCB43fX/qpt45N4+fFsEfKABhpbz7a0rrl9zi/nto9RiEiu1X/hoXu44pTxUHrZ3PGuP4cP6Jn1P0ppejYsSNg/IzMCUcAACAASURBVE9oaGhoydD3ve99r+W84cOHs2VLek5eQAy/IISR6YmBrHv6dkY/XivwRGVHyxTWv7eNxX9fS+MhIxLiwO56Xl1yHN8f9WtKvpxhb6hbBHYRCDZ+sRbKsd2bDwj8inuHjs9vAR9ruF7w+CMJD6PNaugBu89m8d/XAsRk/Juamhg6dCiff/45119/PSNGjAg53tDQwN/+9jdmzJjhq9+lS5cyePDglucvvvgi/fv3j3qcbojhFwQL8Y48SDZOe/oqKwutdYgif9b1V8dssO22BaKtHZAqol0Zvjt3Y4thMWk81MzClSdy5W9X21/kRSAXJ+PXgtfsgXZiu1euaxXr2V13wYPOxttpkmNOaqLNaugBp8/m3bkbYzL82dnZfPTRR+zdu5cf/OAHrF69mkGDBrUcv+666xg5ciRlZWW++hVXvyCkkEyPPHDa09da88s5hlbBb35/PySy70QQy8rwwG570aRTO9Bq1IJXyrlFxs+Gg8bPHJv4/VixGlkzjj7YyNqJ7ZobWsdZ/bXh2n/z1tbKfZGMd6QVfYLK/kb12figS5cujBo1ivnz57cY/nvvvZcdO3bw5z//OS73SBSSwEcQLGR6YiCnvfvgdrf8/rGSyL4TgdvKMBIdu+b7am9J7PPST0PD3RoOthp9MAytXXKfWFLfekkg5FVUV7vbEOa9dpPx3CklcemlRn6A4t6AMn46pRKOM74/Gw/s2LGDvXv3AlBbW8tbb73FgAEDAHjyySdZsGABzz33HFlpXpZaVvyCYEO6Rx644aUcqZ/8/n5JZN+JIJaVYZ9B3Vj99taQtpy8LM4YZ7M3GymW3oo1rK3yecPt3hzwRJluePBmSN+8NXLonK8CPhpWPmVk/XO7f4JW9JE4Y1z/EE8OuHw2HqmqquLKK6+kqamJ5uZmLr30Ui666CIArr32Wo455hjOOOMMAMaPH88993jfrrHu8d91111ccsklUY/VDTH8QpvGKR6/LcXpW/EiuPOT398view7EXTsmm9r5COtDNe/t421y7eFtQ84vaf9FoGXWHorwSvwN29tNfomzQ1GeyTDWvm8c0x99dcwNaBfySsy8vNb7+OITtuYe/MziKeqv7S0lFWrVtkea7QrfuXAkiVLQp6PGjUqqeJhMfxCm8UpHv+rr77i448/bjNx+nZEisP34hXwgp2IL159Jwu7lSEYq3k37LYIADav3sVZdhdEE5serOx3MtxO7cF4zYt/6CBkZQeS77iE4wWTpjH3YBj/WAx9W0UMv9BmcYrH/+CDD8JC3DItTt8Pbgr7WJT3TiK+0ZNvYPTkGzJG1X/8iJ5Ubdwb5rJf/fZWVr+91XGl6HuLwJcbnfgq+/0Y5+YmY+Vv5u2fGiGaJU4x90LyEMMvtFmcXGd2ce1u52cykRT2sRhjNxHf5JlPp62ht2Pz6l2Ox5xU/r63CLyE8QFaQ1NDNof6XEWHYBe6UwrcQg/FXvxOOqq/bi3965Z6FwwvQXCZYCHtEcMvtFmc4vHdzm9rOBnnN/8wnTdm/p5O3Q6n35DT2LRqhe/VeaaJ+NyIJOSzi/92E4/Z5wWwhL05uNGVMrbYv3p0ASXdvkvx2LHGgQsehLnXh6bAzc5rzQjohsdJRwhmvYBImBEIkDrjX7Mb9lcZ7012nlHbPoXV79IdMfxCm8UuHt+JTIrT94NjTH+zYaz279zBx/96o/V8HzH3BR07Urd/v217puG0eg/GetxJPAa45AWwKNwd3Oi5HZrQdXXUPHs3xZ/f3hoff/QZ8MXbtEQEZOeGX+yWec9r1kBoVfxHqrIXfG4qDH/NbuM16cAErOlQ62sU429LWhl+pdT5wAwgG3hSaz3Ncnwk8ChQCkzQWr+Q/FEKmYK5X//yyy/buvfNHPZtTdXvJU+/G6a73s3wr1m6mPqaGttj9TU1rFm6OKNc/U4Cv2A6ds23Xclf+dvvhpw3+45l3jPGObnRNQz4UUBzYDqtqr8ON9qHDhpeADCMbqQyt6WXRi6sE4w5efByfvD2QDLZX9Vq9E10s9Euht+WtMkyoJTKBmYCFwAnApcppU60nPYVcBXwj+SOTshUSktLHQ2f1pqpU6dy4403timjv3DWE0Y4ndYtK3u/7N+5g1nXX82apYttjy+d8wy6qcn2mG5qSttkPU4cP6InZ08c4Lg/n5OXRZ9B3Vj897UtK39zJb/+vdCQvnhkjFNZhss/UP/FnaZD8PK1rSt9u1j9l37amvCn/B5DOBh6R/u+TY+BV+ySDiUau+p/bu0xUFdXx/DhwznllFM46aST+NWvftVybOLEiZxwwgkMGjSISZMm2XoalyxZ0hL3D/Dmm28ybNgwBg4cyIABA7j55psBmDp1KkopPv/885Zzp0+fjlKKlStXxvw60sbwA8OBz7XWm7TWh4A5wLjgE7TWm7XWlUB0/82EdonT3n172dOPFtPtb2f8I+3jW4+vWbqYWddfzSMTxrpOKFLJ8SN6csa4/uTkhf9b7Nm3M5++s9VThj9fGeO8uNG9oJsM4+62Mq/+2si299JPjZTAhV1pyabXdyRhxj84qiAr29s4TJd/MsnO89ceA2ZZ3o8//piPPvqI+fPns3z5csAw/GvXruWTTz6htraWJ5980rWv1atXc8MNN/Dss8+yZs0aVq9eTb9+/VqOn3zyycyZM6fl+QsvvMCJJ1rXwtGRTob/KCD4W7sl0CYIMVFeXk5ubuheaHvb048Wp1S7kZLxBB+3eiHcJhSpxik2f8u6vWHeZBPrSt5u8uCYMS7poXAB71ftbmishfGzDOO+5X1CMwkqOOVyw21fcZ8R4ueVZMf1dyoxXCTBqCzoVEL1vHlsOKecNQNPZMM55VQH8nVES6SyvEoplFKeyvI+9NBD3HnnnS0pf3Nycrjuuutajl988cXMnTsXgE2bNlFcXEz37t1jGr9JOhl+O1+Tv81JsyOlJiulViqlVu7YEZ5BTGhfZHrufT8kIjue3WSibMIV5OQ5ucVDk/VkUu7+aAq45BeFroat2wYdu+Zz9kTjn/vsO5Yx89pFzL5jmbFFYOt2TxINtcYWwUs/tVH7a9iw0PjVryFP9mSmQ1fDa2Gu8LPzoLg31RXLqLr7Hhq3bgWtady6laq774nZ+Dc1NTF48GB69OjBeeed51iW9/zzz3ftZ/Xq1QwdOtTxeOfOnenduzerV6/mueee40c/+lFM4w4mncR9W4DeQc97AVsdznVFaz0LmAUwbNiwqCYPQmZjl5L3xhtvTPWwEo5d1jw7VHY2+R062KryreQXhav0QxIA7dxhlPxtbqbT4d3DwgEzKexPZYXrxCJxqKaJ9e9tCxHuWTPGOVYAnDiS48c+Fhrid+iAZX9aAdowbseNhk9f9patzwvaZSVvGnwncV9hV8Nr4FZ5L1l06Bom5Ns+/VF0XV1Im66rY/v0R1tDJKMgUWV57ZgwYQJz5sxhwYIFVFRU8PTTT8fcJ6SX4V8BHKeU6gt8A0wALk/tkIRMxClVL2ROSt5o69lbDbITSilOOL2MT/9dEXmS4KD78poAKJNy9/s1+mAk3IlU4921NvxvLSF+TuF4ZnvtHsPo1lW7G+5YMVfudjkAcgtb8wdYxwqByAGbcMIk0lhV5avdL7GW5T3ppJP44IMPOOWUUxzPGTt2LLfccgvDhg2jc+fOcRk3pJHh11o3KqVuABZghPM9pbX+VCl1H7BSa/2qUuo04GXgMGCsUuperfVJKRy2kIY4pep1SskbbcGeg6u2s2/BZpr21pPdJZ/OY/pQNKRHzOP3U8/eaYIwsOzssH6CaW5sZNOqFSGpdXGIfqg7cCDq1+E0AUnX3P1e4vnt8JsDwLXdrpqdNUwvXit+J4JX7qWWxENWY26dtLiFEyaRnJISw81v0x4tO3bsIDc3ly5durSU5b311luB1rK8FRUVnsry3nLLLYwfP54zzzyT448/nubmZh599FFuuummlnMKCwt58MEHOf7446Mesx1pY/gBtNZvAG9Y2u4J+n0FxhaAIDjilK3Prj2Sd8DJuB9ctZ29L21ANxiruKa99ex9aQNAzMbfbU882PB7Tcf7xhOP2N5n/66dIefNuv7quK3M3SYddtsB6YKXeH47VBZh7v5goq0A2EI0lf2ipbCrsZoPNtReS+u6hRNW3JfU1X+PG39B1d33hLj7VUEBPW78RdR9xrMsb2lpKY8++iiXXXYZNTU1KKW48MILw86bMGFC1ON1Iq0MvyDEA6dUvXbhe27egf5NPR2N+74Fm1vaTXRDM/sWbI7Z8HvdE/cyQRhYdrbjqttq0N2q6nndenBb5YNh9CfP9L9PWT1vHtunP0pjVRU5JSX0uPEXMe3T2mEm52k81Nyy119QlINGU3/Q3aWum7HN528Sc234ZCrl84qiN85u40zy6t/8fsTzexNrWd5Ro0YxatSolucXXXRRSFy/ydSpU22vt5bzjRYx/EKbwy5Vr1P4npt3wM24N+21d906tfvB65641wmC1zK5ThX7AFvPwjfr1oTk+O835LSImoFoBH3V8+aFrNxMdTYQN+NvFd/pZsMwl116fIgh//c/1oZV8TNpPNTMv57+jHfnbgyr5hdzbXi/RXZiIZZJRqRxJjm1b/HYsXGfILYFxPALbQ5zf97Lvr2bd6DpW2fjnt0l39bIZ3fx6Lp1wauh9jpBiFSCN9Jq/omfXGbrWbDm+A9+7kQ02waJUmcH4yq+CzLOZ10+gLMuH8DMaxc59uVUzS+m2vB+i+xk50Fex0CCIJ+BTbGE43kZZ7Lj/IUwxPALbZLS0tIwQ28n4nPzDmS/Uedo3DuP6ROyDQCgcrPoPKZPzGOPZKhNvE4QzD6dXPNuOoE1SxdTfyByyJ9X+g05zfc1iVZng/80u5HC/hxz85sEKfcbCkp498BEPtn5HWdPgJ3A7rjRRqy9GQYIhqG3iu/85ObPyo0tHM9LMaCkJy0SrIjhF9oFTiK+sWPHMnbsWFvvwMGm7Y7G3dzHT4SqH7yFyrlNELzuyUfSCcQ7yc6mVSt8X5MIdbYVJ0NuTQhn4iXsz1Hpb1G+59Zt5Yysx6graGDD7rOctQJeBXZWnFbhuUXQVBealc9TcYAImOO0KvwhdXH+Qghi+IV2gZuIz6lITyTjXjSkR9wMfbTYTRD8hANG0gnEO8lONP0lQp1txcmQO7V7Cfs7+sAqNpxzf7iwzEb5nptVz+kd/86GurMiewv84hSOZ7cqNwv+BF8X7/umonSvEIIYfqFd4CfEL5hUGfdoE/iA93BAiKwTcDqeW1CAbtZh2wwnnVVuCP6cVP1R7PEnQp1txUu4XXBJ3oKiHFS2c/6ckp0rOXb9P1renxBBosMed6fs1kmR51wCTsl+rNh5C16abN+nboqf+j5aL4WQUNIpV78gJIxMqtAXa1EbPyly7XLumzqBt578g2PynfOuuZ7Rk2+g0+HdQSk6Hd6d0ZNv4NxrrmPyzKf53g2/dOw3GorHjuW4RRUMXPMZxy2qiLtS262wTvW8eaz57igarzybk1+/hSO+fZ+6g40oVEue/vyibAqKjHVUx675DKx6AyyTL1OQ6LTHvb+pdVLkKb7fdKVXfw3o1nA5r2Vx3fbaU1FlL0PYu3cvl1xyCQMGDGDgwIG8++67IccffvhhlFLs3Bn+97ZkyRKKi4sZMmQIAwcO5N577w1pHzx4cMvjrbfeSthrkBW/0C7wE+KXavys2O0o6NjRNge/Uoo1SxeH9OGkE/hm3RpblX5uQQHnXXN9SJ4AO7wKFBNF8OrcKphzO2ZtP2L7CmOlXleHAgrr9zBg3T8A+PaI4eTm53DNI2eF3bvxpW9tq441VlVB+dSwve+G5nyWH5gI+Ijvd0qW8/K1xmo+kms9kgJf1Pe2TJkyhfPPP58XXniBQ4cOUVNT03Ls66+/5l//+hdHH3204/VlZWW89tprHDx4kMGDB7fE8ZvtyUAMfxpQtW0umzY+TF19FQX5JfTrfzMlPcelelhtCj8hfqkmlqI2a5Yupj7oH1EwurnZdq/fTifw5h+m2/bRUFfXIvjzIj5MRXY+x4I4AZyOWY3/u3M3Mvzfj4AllDC7uYF+m17l2yOGc2B3PTOvXdQyUTD7Py3/MArr94SNLaekJGzv21T1b6hzUfXb4WSYzf2HSAlzzLaXr7Xfs2gD6vvXN73OjA9nsO3gNnoW9WTKqVO4sF94djyv7Nu3j7fffpu//vWvAOTl5ZGXl9dy/MYbb+Shhx5i3LjI/7+LiooYOnQoGzdupEeP5G4niuFPMVXb5rJ27Z00Nxuz7rr6raxdeydAWhr/TJ6k2IX4pSOxFLVZOucZdJNzljmvngPd7CxbdxMLxotYNA5OMfn/evozW/W+KaaD8EmB3mW/ci+wGPUDu+v519OftTzf1O/7DFj3D7KbWz1MIYLEoL3vXGBk4OGLwsMi5+yPlDDHbPervveqLUghr296nan/mUpdkzFxqzpYxdT/TAWI2vhv2rSJ7t27c/XVV/Pxxx8zdOhQZsyYQVFREa+++ipHHXWUa9GdYHbt2sXy5cu5++672bFjB0uXLmXw4MEtx1988UX69/eY2dEnYvhTzKaND7cYfZPm5lo2bXw47Qxqpk1S4kGiCvG44RSf32/IaUY+fRdj6MUrEA+lvpcJRLTG209Ugh1uwjgnlb65wrdOGOocVu51+Ye5juHbI4YD0G/TqxTU76Eu/zC+GfRDsnqcRtJVJZFc9n7V92lUiMeNGR/OaDH6JnVNdcz4cEbUhr+xsZEPP/yQxx9/nBEjRjBlyhSmTZvG7bffzv3338/ChQsj9rF06VKGDBlCVlYWt912GyeddBJLliwRV397oq7ePgmJU3sqyaRJSjxIZCEeN+z2x63pcJ2MoZO3IBir58DOQHvBbQIRi/GOVeMQTYU9p2vsVu5NWbls6vf9iH1+e8TwlgmASZ+Xf4F+cyGKJlDZMPQquOj3vsbaQm34hMQWLy57P+p7J21BElPxemHbwW2+2r3Qq1cvevXqxYgRIwC45JJLmDZtGhs3buSLL75oWe1v2bKFU089lffff5+ePUO3bZJp4J0QVX+KKci3T0Li1J5KMmmSEg/ccvUnmoFlZzN55tP8cs48Js98mk2rVjgaw2DKJlxBTk4eTliV9U4RBF5w23pwMt5v/mF6xOiEWDQOYK/Qd8MU09kp6b89YjhrT7ic2vzD0EBt/mF8ceqVHBxY5rl/k7JOf2ZQ/puG0QdjX33lX+C1m9wvdMKLQffisp8+CKZ2MX56iQhw8iCkmRiwZ5G9TsKp3VOfPXvSu3dv1q1bBxiaoRNPPJGTTz6Z7du3s3nzZjZv3kyvXr348MMPw4x+uiCGP8X0638zWVmFIW1ZWYX0639zikbkTCZNUuJBIgvx+MWrMTy640kMO/x8OmR3BiAvq5BcVQDQEnIXvGp2MtBeMrj1G3Iaa5YuZtb1V/PIhLHMuv7qFqPuNF5TYOhm/J0mFF5zABw/oidnTxzgGhJnZuTr2DWfsycO4PgRPR0nDN8eMZx3z/gNi0fN5N0zfsNXHYf4nlwADOqw0P5t/eCvvvppofwew7AHk51nlNVFQXFvGPtYZJe933BApwlHmokBp5w6hYLsgpC2guwCppw6JaZ+H3/8cSZOnEhpaSkfffQRd9xxR0z9mZh7/ObjhRdeiEu/doirP8WYLvJMEMz1639zyB4/pO8kJR7EUognkjbAr3bAyYXfsVNXqqa939JPc30jxxQO5JijB4aNueS24WHXO66itUZlZbmK/NYtX+q4/eC25RDJbe+nBoETZkEcq8Lf6Curxdhbr4FWVb9TGl9zQpGdq2g8FH48Jy8rTCsAoHBKD+he8teRWDPjReuytwsDTMNUvOY+fjxV/QCDBw9m5cqVruds3rzZtt1alje4PVIysXgihj8NKOk5Li0NvUmwkj8nu5isnAIaG/emfJJiV3THTrUfrUDPrhAPGCv+qmnvO/YTSRsQjXbA1hjm5DGow3dbJidunginY44RBId3p2zCFSyaPcs2JwBg224adbvxBuPmto9nDgC/5XCDK+g5TRr6DOrGW898ZmuvO3bNp8+gbnz6ztawSYMmy974q2zfr6uFWDLjReuyz6BUvBf2uzBmQ98WEcMvhGAN1+va7Wy2bXupZZXf2LSXrKxCTjzxkZROVpyK7gAhxj8WgZ41V38wbv24aQOKhvSIeNwOO2M4qOhMjs463vU1mDh5KcomXMH8P82gubGxpS0rJ6fF0A4sO9uIJIggGAxm/66dLeN98w/Tbb0Gkdz2seYAcEvS4xXz/LefX0f9QcPK5+RmsWZ5la3Rzy/K5oxx/Vn897W2noLVNaM5ucP8EHe/1qCGXeVrXHGjuJd9Bb14iwGFtEP2+IUWzHC9uvqtgKaufitbt/7DUcmfStyK7gQTq0CvaEgPSm4bbms4nfqJpA2IVjtgFfx5NfqRygVrrV2fO6X1ze/YybY/06gPLDubC667Ma6pe71grtRNlb6ZpGf9e+Fq7vXvbWP2HcuYee0iZt+xzPacQ3WtVr7uYCNNh+zr29cfbLINCTRZuv9/+KTmfJp1FlobRh8wSut6TbMbT+w0Amnoshfijxj+FFG1bS7LlpVRsehYli0ro2rb3FQPyTZcD+z/yaVaye9WdKeysrLlebwEen76cVpdm+2RjnvF6fysDjkh9+oy/jhHT4Jdwh/d1BQSLTCw7GxGT76Bgk6thj47L48effrZ9tlvyGlh11pz+icyo59TAh8zSY+JlwnC28+v87UFHymMcOn+/+Gt6p/TqPNRKqCh9JtjP16UXmqI/4p740kMKLQZxNWfAtI1EY4fY55qJX9xcbGj8Q92+cci0LOe77UfO21A8Ko70nGvOPVTPLa/5zwDfkLnGutblWz1B/bz9eqPba/dtGpFyPNkp+51Mr7WdrcJgunmN1388eT0jn8nN8syxlTFwYvLvl0iK/4U4JYIJ5U4G/PQGCQ/Sv5EeTbKy8vJzc21PRbs8u88pg8qN/RrHq2R9dpP0ZAedBl/nOOqO9Jxr8SjH6+hc3Zhf07EIzNgLDiF8VnbvU4Q4k1w+d0Q0iwOXmi7yIo/BaRrIhyncL2ePceze9di3+GGifRsmAK+l156yfa46Q2wCvSiTbvrt5+iIT0cFf/BfRz2oxNiygLodJ9ItGTrcyi7a92D92PMvcbbJwpTYGdV41sr3jll6wueIBQU5VB3sDHsnGgoKMqh7NLjUf+JQVQnpJy9e/dyzTXXsHr1apRSPPXUU5xxxhktxx9++GFuueUWduzYweGHt/4tLFiwgFtvvRWAzz//nKOOOorCwkJKS0uZNGkS48aNo2/fviH9nHvuuQl5DWL4U0BBfklAQBfenkrinVMg0Sl+S0tLW8L5rBQXt2ZEj9Y4WrHrx0+oYDxSAMejdoA1nW4wZhif1TXvJRUwJF645wWvIXxeJghllx5PxTNraG4K17rk5GWRk5sVcWIQdv/CzIiDb7PU7Ib9VdB0yEh41KkEOnT1fHm0ZXnHjBnDmDFjACNu/+GHH2bYsGEA6ZurXynlVGBYA3Vaa+/xPu2cTE6E46c6n7NnI3zSEy3l5eW88sorNAeFjGVlZVFeXh63ezjh15BHE8YXy/2ccHLbdzq8O5NnPh3SFskzcNJZ5WxatSLmePt4ExyP73YOuE8QnJL6WMvwOin5O3bN58rffje0MYPi4DMeaxXBspvgqKGtmZmaDrV6XzwY/3iW5U0lflb8m3GSeANKqX3A08D/01rHxzfWRknXbH2RXPN+XfdOng1QVG2bG7fXqyx5UK3PveA1GVAwToZ8z/Pr2PPPdWEr8lgjDGKdOJh4FfRF4xnINLxOEJzOWf/eNnJy7TP1gYteQER1iceuiuD826DsZjjuvNbzdLPhAfBg+ONZltdKMsvy+hH3XQZsAe4Czgs87gK+AiYBU4H/Bu6O7xDbJiU9x/Hd7y6l/JzP+e53l6bc6ENk0aFfUaLhwbAzwjpuQsaKigqaLOFoTU1NYfH8bpjJgMwtAzMZUHBYoB2OBlu3Ht/70gYOrtoOxB7GF6/QRK+CvkWzZ7l6BjLd6MeKGQ4YydVvlxtASAJ2KYkb62HFk+HnNtnkXrbBLMv7s5/9jFWrVlFUVMS0adOoqanh/vvv57777ot6uGVlZXz00Uctj0QZffBn+H8G3Ki1fkBrvSjweAD4JTBJaz0D+DnGBEGII8mK+Y8kOvQrSjQmM4nNA+AWz+8Vr8mArHgx2LqhmT3/XMfW+94lf8BhMUUYxCv+3ykhj7k3v2bpYp74yWWOqXpTrdpPF9yS9QSz9Pn1SRiNEIZTlMSB7eFt2c4VLYOxK8v74YcfhpTl7dOnT0tZ3m3b0nPS58fVPwL4xKZ9NWBm7HgXiFqaqpQ6H5gBZANPaq2nWY7nA88AQ4FdwI+01pujvV8mEA9lvNu+fPAxYx4YHrdsig6jESUW5B+ZUCGjUzx/sLgvEtFOHpxy+dvRXNNIzfJt5PbvTPOueltxXiThXrzi/91y4bu5901SrdpPF7yG/cUrKkDwiVNK4o6WbTGVZQj8PBBclveEE04IK8tr0qdPH1auXBmi6k8n/Bj+L4HJwC2W9p9iuPsBugO7oxmIUiobmImxhbAFWKGUelVr/VnQaT8B9mitj1VKTQAeBH4Uzf0yhViV8XYTh88+u5X1639NY+Mey9nhRj9YdBiNKDHRQsby8vKQnP0Aubm5vsR90U4eTKO855/rPN+rYeM+2xA+L8K9eIUmgnNSHS/x+qlW7acLTuGAQprgVEVw1G3GCj9KVb9ZlvfQoUP069ePp59+OvJFHrDu8d91111ccsklcenbih/D/0vgRaXU94AVGD7c04D+wA8D55wGRJt3cjjwudZ6E4BSag4wDgg2/OMwtAQALwBPKKWUtiYXb0PEfER0AAAAIABJREFUGvNvn4a3wcboB2Psy1u9A9GIEhMtZDQFeH6FecHEMnkwi+742We3E+N5Fe7FKzQxmBbl/q6dQQnk7cnv2Ml1bz+4r3RS+SeCPoO6sfrtyBEqOXn+xaZCHEhQ9EQsZXlNlixZEvI8bcvyaq1fV0odB1wHnIBhHV4F/qS1/ipwzh9iGMtRQLBfZgvG9oLtOVrrRqVUNdANCNl0VEpNxvBO2MZTZhKxxvxHt5euKT9no+2RaEoIJ7rscGlpqS9Db3c9RD958OPyB3sxXryEe37x4to3ycnLp/yqyZ772r9zBwtnPQHQJo3/5tW7PJ3X1KRZ/94239UBhTgg0RO2+Ergo7X+Grg9QWOxl3/7Pwet9SxgFsCwYcMy2hsQq6vcOaROCCaWyYPVBZ/VIYfmGud9XTsxXrxqCvjFayregk6dOOfKya4G3K6vxkP1LJ3zTJs0/F7d/LqJkPz/gpBqfBl+pVQHYDDQA0tEgNbaPneqd7YAvYOe9wKsFss8Z4tSKgcoJkpNQaYQq6vcbuIQiZycw6Iaa3vG6oI/uGo7e15aDw2h804nMV68hHt+Xe2uCn2lfLnr/RT8yXT8huiJFkBIJ/xk7jsXeA7DtW5FYyjxY2EFcJxSqi/wDTABuNxyzqvAlRjRA5cAi9ry/r5JLK5y87r16+6jsWlvxPOVyuX44+OfisFPxr9MxqrMzx9wGPVr90QU48VDuBeNq90pFa9dFr9IOPbVBqMArCV+I+FUOEgQUoGfOP4ZwOtAL611luURq9EnkO3vBmABsAZ4Xmv9qVLqPqXU9wOn/QXoppT6HLgJuC3W+7YHSnqOIyeng8PRLHKyuwCKgvwjGTjwwbgbZDOywNhy0C0hiYnKR5AqTGW+6bJv2ltPzfJtNNc3ctiPTqDktuGuhrxoSA9KbhtOr2llEc+1w83V7kSkmH4/xLOvdMdtBZ+TlxX23FogSBBSiR9Xfx/g+1rrhG0Ya63fAN6wtN0T9Hsd8F+Jun9bxlnkpznrrA8Seu9EF+tJF+yU+QC6timqnPp+icbV7hbT75d49pXuuFX2O2Nc/4gFggQhlfgx/Msw1Pz+fFxCSrC61nOyi21d/cmoCJiuZYjjjZsCP5qc+l4xtxc6ZHWipmlf2PFO3Q533ft3iumPhnj2lc64Vfbzkv9fyFyiLcsLRhifU/nd7OxsTj75ZBobGxk4cCCzZ8+mQ4cOdOzYkQMHDsT1Nfgx/H8CHlZKHYmRwS8kx6nW+sN4DkyIHrukPUrlArkEf2zJqgiYrmWI442TMt/E7Vi05XaDE/+UHjaSFbvm0xRUIysnL59+Q05rV2F2ycCtst/697bJir8NE21ZXhOn8ruFhYV89NFHAEycOJE//elP3HTTTfF/Afgz/C8Efs6yORYPcZ8QJ+xc61o3kJNzGDnZhUkX2GVyGWI/RIrndwrNi6XcbvD2wjGdTgKgcs/b1DTto0NuMaWdy6hcspTGhviH2bWnZD122K3szcI9pifgwO56Fv99bcv5QnJ5ZdU3/G7BOrbureXILoXcMuYELh5yVNT9Jassb1lZWcQiYbHgx/D3jXyKkA44udAbG/dy1kj3jFOJIF3LEMeboiE9qP+ympr3toVll3ALzYul3K7Vi3BMp5NaJgAmy3eGry4gtjC7tpysJ5YVu13hnsZDzRLHnwJeWfUNt7/0CbUNRiryb/bWcvtLRrmZaI1/PMryRiq/29jYyJtvvsn5558f1Ri94Cdz35cJG4UQV9LRtZ7o7H3pwMFV26n9YHu40S/Mpsv3j3U04tFm7Tu4aruR0ipCQGuH7M6Oe//R0laT9cS6YndS+0scf/L53YJ1LUbfpLahid8tWBe14TfL8j7++OOMGDGCKVOmMG3aNG6//Xbuv/9+Fi5cGLEPJ1d/bW1ty4SgrKyMn/zkJ1GN0Quuhl8pNR6Yp7VuCPzuSBwS+Ahxor241tMNJ1W/Uop9Czaz55/rbPfvo8nat/uVDdQs95ZExmnvP5Ywu7aarCfWFbub2l9ILlv32ictc2r3gl1Z3mnTpoWU5QVayvK+//779OzpzdMTvMefaCKt+F8AegLbad3jt0P2+NOItuhar6ysjKkQTyKwCvKcVujNNY0QSOFrt3/vlLUvf8BhVE17P0zwd3DVdmejryCrMDRlcMvef/VSahr3edqPj7R/31aT9cS6YndT+wvJ5cguhXxjY+SP7FIYdZ/toiyv1jrL7nch/WlLrvXKysqQ6nnV1dXMmzcPIGXG306Q5xXr/r1d1r78AYdR+8F2W8HfvgWbXTqH4rH9wyYSfbqezOBrxnuKFPCyf1824Yqw4j5tIVlPrCt2N7W/kFxuGXNCyB4/QGFuNreMOSGmfmMty5vM8rtO+MrVL7RvUpV2t6KiIqRkLkBDQwMVFRUpM/xObn2vWCcK1lz/VdPedxT8uU0ysrvkx5z+18v+fVtN1hOPFbvE8acH5j5+PFX9EFtZXrfyu06x+vGO4Qf/RXp6A2XYF+n5fRzHJaQZdrkB1q69EyDhxt/pDyWZ9auteFrhuwjvIlXdcxP8uW0rmJED1omEH7zu37fFZD2yYm9bXDzkqJgNfVvET5GeicBTQCOwg9B/aRoQw9+GSWXa3eLiYlsjX1xcnND7uhEpWQ/gaPS9VN1zKu2b1SHHMV9Ah9N7xiUzYFvdv/eKrNiFto6fffv7gEeAzlrrPlrrvkGPfgkan5AmpCrtbmVlJYcOHQprz83Npby8PKH3dqPzmD6oXP+yl+wu+XQZf1xEA+1UdFJrTdGQHnQZf1yL10AVZpPVIYea5duomva+EeYXA+2p2I4gtEf8uPqPAJ7UWjdFPFNoc6QiN4BV1GdSWFjIBRdckFJVv2m49zy/LmIcfTAltw33dJ6utf8zM9tNV74pMmwOCJj8ZP1zoq3u3wttA601SqlUDyNtiKYyvR/D/wYwAtjk+y5CxpOK3AB2oj4w0mTG0+hHGypoGla3NL3BZHXIsQ3Ps8NrbH8sWf/caIv790LmU1BQwK5du+jWrZsYfwyjv2vXLgoKCnxd5yWBj8m/gAeVUidhX6RHEvi0YVKRGyAZor5YQwWtCvqsDjnG6rvBMgvPVjTXOcfzW3MCWMP5wF4bEG3WP0HIRHr16sWWLVvYsSNcg9JeKSgooFevXr6u8ZLAx8odNm2SwKcdYDX+mzY+HNIeb5Ih6otHqKDV7W41+lkdctBah7nvzZU5EJYToPaD7RQO7UH92j2uHoJosv4JQqaSm5sbUtJWiA7PCXwEIdkhfeXl5WF7/H5EfV5c+PH0Kjim7M3LptllZe7krq9fuyeiJsAp61+kqAFBaG+8vul1Znw4g20Ht9GzqCdTTp3Chf0uTPWwUoIYdsEzbiF9iaC0tJSxY8e2rPCLi4sZO3asp5W46cI3DbjpwreWunTyHkTjVYgUe2+HW1igF3e9VeGPavUkxKruT1eq581jwznlrBl4IhvOKac6sDUjCE68vul1pv5nKlUHq9Boqg5WMfU/U3l90+upHlpK8BPH/xTwqdb6EUv7TcCJWutr4j04Ib1IRUhfaWlpRENvt7L36sKP1asQjJvb3W1l7pSNz6u73k5kGA91fzpSPW8eVXffg66rA6Bx61aq7r4HgOKxY1M5NCGNmfHhDOqa6kLa6prqmPHhjHa56vez4v8esMimfVHgmNDGcQrdS2W5X6eVvVcXfixeBSt2sf2mcbeuzIPj+d2u84qbur8tsX36oy1G30TX1bF9+qMpGlFm8fqm1xn9wmhKZ5cy+oXR7WbFu+2gfVErp/ZoyZT31084XxfALmnwQaBrfIYjpDPpWO7XaWWvlLKNb7Vz4XvxKnghUo58pzS6sebWh/aj7m+ssvcuObULrZjubnPla7q7gTa/6u1Z1JOqg+HfkZ5F8cvQmEnvrx/Dvx5jZT/D0n4h8HncRiSkLelY7tdpZa+1Jjc3Ny4ufD9EmyM/ltz60H7U/TklJTRuDU8klVOSOq9TptCe3d1TTp0SYpQBCrILmHLqlLjdI5PeXz+G/xHgT0qpHrS6/MuBXwDXx3tgQnqSLuV+zX19J5RSnHLKKWzYsMF3Yp5MpL2o+3vc+IuQPX4AVVBAjxt/kcJRZQbJcnenI6bhTaSqP5PeX8+GX2s9WylVANwF3B5o/ga4SWvtryCxIMSAUyrfYLTWfPzxx1Hv12ca8dguyARMAd/26Y/SWFVFTkkJPW78RbsX9nkJVUuGuzudubDfhQldeWfS++urLK/W+s/An5VS3QGltW6b8UJCWuOUyteK30Q8qcKatS9agx3rdkGmUDx2bLs39MF43VtOhru7PZNJ768vw2+itZZ8iULK8JNcJ57pfROBme0v2jC8eE0ahMzF695yMtzd7ZlMen/9xPF3Be7H2NfvgSUUUGvdOb5DEwR7nFL5Op2bzsRSZCfWSYPQNvCzt5xod3d7J1PeXz8r/r8AQ4BZwFZ8FSN1JzCp+CfQB9gMXKq13mNz3nzgdOAdrfVF8bq/kFnYJd3Jzs5Ga01zc6sR9aLij7YyX7yIJQwvUZX5hMwik/aWhfTAj+EvB87TWr+XgHHcBlRoracppW4LPL/V5rzfAR2A/0nAGIQMwTTMVoNt12ZnxIONfTB+K/PFg1jC8NpL7H60VM+b5ygCdDuWacS6txxJGCg57tsefgz/duwT+MSDccCowO+zgSXYGH6tdYVSapS1XWh/OCXd8ZLe1y0iINmCwFjC8NpL7L5Xgo15dnExTQcOQKNRBjk4tS/QptL+xrK3HEkYmElJaQTvKLvsZrYnKvUj4FLgSq11XCcASqm9WusuQc/3aK0Pczh3FHCzm6tfKTUZmAxw9NFHD/3yyy/jOVwhg5k+fbonfcDUqVMTP5gA0Qr0rHv8YEwazFTA7QlrDn8nco48EsA+CdCRR3LcIufcEG0JcxVvt0UAUFJUwsJLFjL6hdG255jHhfRFKfWB1nqY3TE/K/67MPbgtyulvgRClkxaa9clklLqLcBu0+lOH2PwhNZ6FoYWgWHDhsVNiyBkPl6MfrIFgV7D8OwmCF3GHyeqfuxz+Nvhltq3vaT9ta7i7TCFgU7CQbvJgNOWgGwVpB9+DP8LsdxIa32u0zGl1LdKqRKtdZVSqgRjW0EQ4k6kiIBkpPWNBicFf5fxx1Fy2/AUjy71eDXaZmrf9pz21y78z4opDHQSDoJh6E0DbrclcOc7d3Lb0ttCrpGtgvTAc3U+rfW9bo8Yx/EqcGXg9yuBuTH2Jwi2lJeXk5uba3sslsp8iaa9VN/zSvW8eWw4p5w1A09kwznlZHv00uiaGjqeNRJVUBDS3p7S/npJITuy10gAV4HgjA9nhPxunUw06Sbb6+qa6njgvQe8DFVIEH7K8iaSacB5SqkNwHmB5yilhimlnjRPUkotBf4PKFdKbVFKjUnJaIWMxa4M7/jx45k6dSo33nhjWhp9EAV/MOZ+fuPWraA1jVu30nTgAMoyoVO5uVBYGNLWtHcv1S+/QvEPLjb2+5Ui58gjKfn1fRkp7IsGL2F+CzYvANxX5cETCL/56KsPVadtydr2gJ8EPnkY+/GXAUcDIX9lWuvsaAehtd6FES5obV8JXBP0vCzaewiCSbzK8CYTUfC3Yruf39gIXbqQ06FDSIje9umP0lhbG3KqrqvjwL/fznghX7R753bhf1b21u9t+b2kqCRingC3LQEn0rFqXXvBz4r/1xhu+EeAZuAWYCawC7gu/kMTBMGk85g+qNzQP9e2WH3PC077+XqvYayOfOhBjltUQfHYsY7nZqKQ7/VNrzP6hdGUzi7lzOfO5O5ld1N1sAqNbtk797KKvrDfhUz9zlRKirxpGqacOoWC7NCtEWueALtzIpGOVevaC34M/6XAtYFCPU3AXK31z4FfYbjnBUFIEEVDetBl/HEtK/zsLvltPmzPuo9fHUiu5CbCM2PyI52baUI+UzxnGvrqQ9U0NIfmojDz83vhwn4XsvCShRTnOWsjzElE8ERBoSgpKmHqd6aG1QEIPqdDToeIY4hnZsHgSdHoF0bLNkIE/MTx1wADtNZfKaWqgIu01h8opfoCH6drrv5hw4bplStXpnoYgiB4pHrePL69/7c07d0b0q4KCij59X0AEWP2zZh8p/j+7C5dOOLOOzJmX98pnt4OhfLs+n990+thynuTWGP1f7P8N/zf+v+jWTeHHSvILgibPESLXXhiPPvPVNzi+P0Y/rXAVVrr5QGR3Zta698qpS4Hpmutj4jfkOOHGH4h3ZEKe61ESsQTbNC3T3/UNiwv+NyOZ41k/5vzwyYRwZiTACBiGt/1723j3bkbObC7no5d8zljXH+OH+F/5eo3ZXDp7FK0z/IoXo3fybNPtm1XKCqvrPR1Tyei1SN4uU6SDNkTL8P/AHBAa32/UuoS+P/t3XucXGWd5/HPL51O6JCxmwhJOg0hgCFEQAm2chOUm1HZkAiMzgyuYUeWVYdRcIgmL3CGhddMohHBGcZVxnUBRWUFTIJxN0rCRSHARAPhEsMlIiQdLmNIuyFNrr/9o041VdXnVJ1TXZdTVd/369Wvrq5zquo5ufTvPM/ze34PPwI2AT3AYneveCGeSlDglzRT9b18z55xZtFgjhnT1z8d//y42trAbLDEb/azcGfkpEmMv/wyXhn/Xu697Xfs2fXW39XIUSM4/cKjEgX/sJub7GhGVPBP0uPPFSf4pTVwxu3JR90UVfLGpREVC/xJ1vEvcPd/DB7fAbwf+BfgvLQGfZG00/r8fKWS7grn5sdfftmQNfll2bs3P+gDBJ2ibN7Ahn+5LS/oA+zZtY/VS59P9FFhqxL8zTd59fobIl8Tljw30kbSNboLwyJfFyeBLuq9B/YMJJozD5tnH87ce1htgLA8hqhcAe1OGC1W4DezdjO73cyOyD7n7o+4+zfc/WfVa55Ic9P6/HzFku7Ciux0zppF97XXDNbgrxZ/8016nrgz9Nj2rcn+rspZaZBNnpswogtzOLDfufT+Dn62/5dYN3ddZIZ+nOBXmJjXOaoTM2Pbzm2xVwwUJh9ueWMLX3nwK1z166vKWnkA0Tcthc/HWXUg+WIFfnffDXwIEk4yiUhRUevwW3F9PkT34K2rK3IovHPWLKauWln14L/fztfzfp7wyqOctPoqTr/vb/JWHZRS7kqD9z+1jxuv287ti/bwrW/t5eQHXx9cwTDc4JfN8l83dx1j2sckXjEQ1jvfvW83ezx/FCXJyoO4Pfk4qw4kX5LlfHcB51WrISKtSOvz8+X14IOqepMWf42jHl5dMgN/WMP+bW0wsng9M3v7BEaOyvxdTXjlUY7a8EM6dr6OMXQZYdJ2xikZXGyKoJLBL25PO+6xcs9NcjOTe+Pyiwt+oaBfQpJNel4ErjKzU4E1wBu5B939G5VsmEgryCbwKav/LZ2zZpW9zG7EfvuxN8YufYW6Pv7njDn++MiVArbffnTP/ztOH38Uq5c+z+Grl9FW0CvOBuFSbc8ez83qH/uB03j1+hvo+9KXI7P8S00RnHP4ORUJeFFV+IpNGySp3Bd37j17LdrZr/KSZPX/vshhd/fDK9OkylJWv0hji7P0LXQZYE5W/t4dOwYr+4XJLhOM+5nrp79zMPkvT8Gqg7jXFyfLP2oFQ2Hbh6ucdfFhr2kf0Y675w33a3197RTL6o/d43f3w3LecGzw3PbhN09EJFxhUMwOqQN5QTG0fn8Q9IsV8skq7E2XGnUY2d1dsa19iw3h57Zh/OWXhd4gVHpXwXJ62lGvSfo+Uhuxe/wAZnYZ8EUya/cB+oBvADd4kjeqIfX4RRpX3F5unB54/9130zd/QWbpXon3K9XjL2ctfpQkowdJC/9I66pIj9/MvgZcAiwGVgdPnwT8PdANfGmY7RRpKarYV1rcpW/FeuC5wdI6O+GNN/Ddb83PF/aa44wyhM3TlxuEk4weDCf/Ia5yq+xJ40gyx78VuCQo3pP7/AXAd9z97VVo37Cpxy9ppIp98cTt8Uf1wDs/Nof+ny7JH0ofOZK2sWPZ298fGrCfOfGk0BK/lZ5LL9X2ckYPhitsrh6gc1QnC05YkOgGQDcQ9VWRyn2BsPqH68p4H5GWpop98RRb+pa7e9+r199A58fm5C0D7L72Grbf/8DQef09e7AxY5i+/unB7Xuz+u++O7Kuf7W28g1bwpgN+lE7FFZL2Hp8gP5d/YmK74QV9EnyeqmuJD3+G4Lzv1Dw/PVAW7BFb+qoxy9ptGn+ryKPHbzo1Bq2JN0Kd+qzri66gw114vSSk2bfF6v9X60ef5Rioxjb739g2FMMYT3yBb9aUHQzoM5RnYxpH1OyF5/W+v+tpOw5fjP754JzP2lmM4GHg+dOACYBt1WioSKtoq1rdGhZ3lat2BcmLPB5fz99876UKbhTkKQXlgmfNPu+WK++0tnzpURl+2/70Y8Hf45a5VBK4ZB+tkfeObqTbTujlz327+qnf1d/3muAIcG/nCJA2XYNZ3pA0wvxlBqiPzbnazrwG2ALcGjw9TLwW+CoKrYxde58eSu9Dz1F972P0fvQU9z58tZ6N0lq4I21r7Jl0aNsmv8rtix6lDfWvlr2e6liX2lRS/SA0Mx8GBq4k1bJi7ohsK6uyMBareH4uFMLpTb4CRO1AY67D6mWV0xUCd5yNs4Z7vSAphfiKxr43f30mF9n1KrB9Xbny1u5YsNLbNq5Gwc27dzNFRteKhr8daPQ+LLJeNle+t5tO9l217NlB//9Z4yn67ypgz38tq7RSuwrUM6cemHgLjZ/HibqRiE7vVAoOyqxp68P3BOV7s19j7AbhyQ1AZL+WUX1vP+060/MfsfsYb9XOXsHFNuNL84uf3F385NkJXsFWLhxCwP78ufABvY5n1//IgDnTxwHZIL9wo1b2LRzN8ZbuxtlbxRyz5X0K5aMV26w3n/GeAX6IqKG6aNE9eSTLIFLukwvbvGdKMWWDoYV7MlWIyyUtHBQsbK8D2x6IPQ1hoXO/4f14sspAhR1M5LtuRdOS+R+TrHXJ9lHoFUo8Ce0eefu0Of3Apeuf5FH+zPFDG/t2zr4X6Twv8rAPmfhxi0K/A1E2+fWXmjgK9TWBvv2VbSYTZIbhXK22M1V7MYhm0hYWNO/cHliOdX7vnD8F0LL8mYT/MI4mWmAsNeESbp3QNTNyAgbEdmTz33/cvYYaFUK/An1jG5nU0Twd+CWvnjD+FE3EJJOSsarrGceeZnVS59n+9adjB03mpNmH8GRJ+T/gs7rfUdtnFOHte65hlu6t9SNQ9hNyOBmQsPI6i/WI//mb78ZmZH/heO/ULXkuaibkbDlhTC0J1/sZkbyKfAntODwbq7Y8NKQ4f6keka3V6hFUgtvmzkltOCOkvGSe+aRl7n3tt+xZ1fmz3L71p3ce9vvAEKDfzaopbFc7XDr55dz45C0et/yjctZ+MjCwWz8rtFdzH/f/MgeebEAWqkdAMNE3YxE3YgU9uS1m198CvwJZYfnP7/+RcLzikvrGGEsODz5Zh5SP9o+t3JWL31+MOhn7dm1j9VLnx8S+HPVolxtUsMt3VvtjXeWb1zOVb++Km+HvG07t/GVB78CDF2Gl/tcPQJo1I1FWDXBHbt3sHzj8rzzq3lj0kwSbdLTiKpVwOfOl7dy6foXi5S6CHfw6HYWHN6t+X1pWf/6mVWRx/7m2y2zQGhQNUcyogrpQGMV01m+cTmLHl00pMaAtvmNVpFNeiTf+RPH8Wj/9rwkvigGfGrSOL46bXItmiaSamPHjWb71qH5EmPH1T9foh7TCZUcySgsYBMV9KGxst2zuQeFgT8syU9KS0WNfTMbZ2a/NLNng+8HhJxznJmtNrOnzGydmX2iHm3N9dVpk7lx+mQOLjJf3wbcOH2ygr5I4KTZRzByVP6vnpGjRnDS7CPq1KKMSqzJr6ewAjbFhGW7x1kvXy9arlc5qQj8wHxgpbtPBVYGPxfaAXzK3Y8GPgzcYGZdNWxjqPMnjmPNyUdjEcf3ofX6IrmOPGEip1941GAPf+y40Zx+4VFF5/drodjSukqo9oY7ix5dFJkBX6h9RPuQbPe0V74rpxqghEvLUP9s4IPB41uA+4Av557g7s/kPO4zs1eBg4DowtI1FLXMT9n7IkMdecLEugf6QsNdk5+rcMqgcP19uTX2oyzfuLxojf3OUZ2hWf25ilW+S8NQupbrVU5aAv8Ed98C4O5bzKxoqrSZvQ8YBTxfi8YVyq3K10ameM8BbSNoN2N3TrKksvdF0idqHn+4a/Jz37+wGl/uxjpZSSr8lVKsLG3cJL60D6VruV7l1Czwm9k9QNgt/pUJ36cb+D4w1933RZxzCXAJwOTJlZ1bz9bqz67jzy7pe33vPtqBA0a2sW3PXnqUvS+SOklL5JaztC50c6EI5YwmhCkWnOP2iNNU+S5ql72w5XrakS+5mgV+dz8r6piZvWJm3UFvvxsI3fnEzN4GLAeucveHw84JPusm4CbILOcbXsvzhdXqz9oN7N82gvWnHlvJjxSRCklaIrecrP4kwTzpaEKUqKDdOaozdhBMy1B61JbBMLTuQJJz5S1pSe5bBswNHs8FlhaeYGajgJ8Ct7r7T2rYtjylSu2qFK9IesUpkTt11Uqmr3+aqatWljUMHxnMLT8FuJzRhKis+6jd8BacEF53P8w5h5/D1SdfTff+3RhG9/7ddVkjn2SXPe3IV560zPEvAv63mX0aeBH4cwAz6wU+4+4XAx8HTgPebmYXBa+7yN0fq2VDi9Xqzx4XkaHi1OcfrlLr8Cs1j19M1JRB58fmsP3+B8oeTYjTux3ukHc1K9/FHZJPkmuQ9ryEtEpF4Hf3PwJnhjy/Brg4ePwD4Ac1btoQxWr1K5lPJFyS+vzlKjZ/nw2w1S6Rm/tZlS4EVCrrPs3lapMMySfJNUhTXkIjSctQf8M4f+I4vj7tkMFmatpcAAAUK0lEQVSiPW3B8wePbufr0w5RMp9IiGL1+ZMothY+zjr8zlmz6L72GkZOmgRmjJw0qSo7/FViyqBQI/dukwzJR01bhOUaJDlX3pKKHn+jOX/iOAV4kQTCSvQWez5MqR593HX4adzsJ45G7t0muWlJMm2hJX7lUeAXkaqrRH3+Yj36Sq7DT6u0ZN2XI+lNS5JpizRPcaSVhvpFpOoqUZ+/VI9+/OWXYfvlD/tWev6+ntKSdV8ODcmni3r8IlJV/XffjV1/A6f1bWFnxwE8N2UWb0w/NXFWf6kefbWS6tIkbu82bUVtNCSfLuZe0fo2qdPb2+tr1qypdzNEWlLhvDxkeuHlJNRV8r2aWWEGPWjf+lZkZr9x996wYxrqF5GqqeSOd7XKyG90KmojpWioX0SqppI73kF6MvJLFQqqp0Ze9ie1ocAvIlXTjJn2cQoFLVm7mcUrNtC3bYBJXR3MmzmNOTN6atK+Rl72J7WhoX4RqZpKZdoXK9xTa6WmL5as3cyCu55g87YBHNi8bYAFdz3BkrWba9I+ZdBLKerxi0jVVCLTPk4Pu5ZKTV8sXrGBgd17844N7N7L4hUbatLrVwa9lKLALyJVNdx5+VKFe2qt1PRF37aB0NdFPV8NKmojxSjwi0hqhO3gt7fCCYIwvOS8Uhv9TOrqYHNIkJ/U1VF2e0UqSXP8IpIK2R38sqV9B3fwGzc+9PxyEwSzUwd7+vrAfXDqIG7eQKllhfNmTqOjvS3vNR3tbcybOa2s9opUmnr8IpIKUTv4PXfYLKa+8YOKbaVbiamDYtMX2Xn8emX1i5SiwC8iqRC1U9+LY2dw2rVHVmzdfKVrC4SZM6NHgV5SS4FfRFKh2A5+nbPOqFgiXzPWFhBJQnP8IpIKldjBL45m38VPpBT1+EUkFbI79RVm9SfZwS+OVtjFT6QY7c4nIiLSZLQ7n4iIiAAK/CIiIi1FgV9ERKSFKPCLiIi0EAV+ERGRFqLALyIi0kIU+EVERFpIKgK/mY0zs1+a2bPB9wNCzjnUzH5jZo+Z2VNm9pl6tFVERKSRpSLwA/OBle4+FVgZ/FxoC3Cyux8HnADMN7NJNWyjiIhIw0tL4J8N3BI8vgWYU3iCu+9y9+wOHqNJT9tFREQaRlpq9U9w9y0A7r7FzMaHnWRmhwDLgXcA89x96BZbIiJSMUvWbmbxig30bRtgUlcH82ZO05bDDa5mgd/M7gHCdtu4Mu57uPtLwLuCIf4lZnaHu78S8lmXAJcATJ48ucwWi4gMT6MHzSVrN7PgricY2L0XgM3bBlhw1xMADXUdkq9mw+Xufpa7HxPytRR4xcy6AYLvr5Z4rz7gKeDUiOM3uXuvu/cedNBBlb4UEZGSskFz87YBnLeC5pK1m+vdtNgWr9gwGPSzBnbvZfGKDXVqkVRCWubJlwFzg8dzgaWFJ5jZwWbWETw+ADgF0L8+EUmlZgiafdsGEj0vjSEtgX8RcLaZPQucHfyMmfWa2XeDc6YDj5jZ48D9wNfd/Ym6tFZEpIRmCJqTujoSPS+NIRXJfe7+R+DMkOfXABcHj38JvKvGTRMRKcukrg42hwT5Rgqa82ZOy5vjB+hob2PezGl1bJUMV1p6/CIiTWXezGl0tLflPZckaC5Zu5lTFq3isPnLOWXRqrrkBsyZ0cPC846lp6sDA3q6Olh43rFK7Gtwqejxi4g0m2xwLCerP03Z9HNm9CjQNxkFfhGRKik3aBZLDFQQluHSUL+ISMo0Q2KgpJd6/CIiNVaqsE8zJAZKeqnHLyJSQ3EK+ww3MVCkGAV+EZEailPYR9n0Uk0a6hcRqaG48/fKppdqUeAXEamiwvn8rjHtvL5j95DzoubvG32jH0kfBX4RkSoJW48fJmr+Pk3r+aV5KPCLiFRIYe98x649Q+bzCx0wpp1/mHV0aCDXen6pBgV+EZEKiNu7LzRm1MjIIB6VD7B52wCnLFqVqBKgpgskS1n9IiIVENY7j6NYUZ5i6/bDlgGGibN8UFqLAr+ISAWUW1WvWHAPW8+fq3AZYJhi0wVp2AhIak+BX0SkAqICeFdHOz3BMSs4VqooT+56/iilbjiKTRdoJKA1KfCLiFRAVLW9q889mgfnn8ELi87h+k8cl7goz5wZPTw4/4zI4F+qjG/UcTNKFhKS5qTkPhGRCoizDe9wivLMmzktL3kQ4pXxDXtd+whj9z4PPV8bATU/BX4RkQqpZrW9ODcWcV+3Y9ee0CJCoI2AWoG5h9/1NYve3l5fs2ZNvZshIpIah81fTtRv/hs+cZyW+jUBM/uNu/eGHdMcv4hIiymWiKig3/wU+EVEWkyxRERpfprjFxEpIk7Vu0arjFduvoA0BwV+EZEIcTbJadSNdLTtb+vSUL+ISIRiVe+SnCOSJgr8IiIRota05z4f5xyRNFHgFxGJEJX9nvt8nHNE0kSBX0QkQlT2e261vDjniKSJkvtERCLELcNb6hyRNElF5T4zGwfcDkwBXgA+7u6vR5z7NmA98FN3v7TUe6tyn4iItJpGqNw3H1jp7lOBlcHPUa4F7q9Jq0RERJpMWgL/bOCW4PEtwJywk8zsPcAE4Bc1apeIiEhTSUvgn+DuWwCC7+MLTzCzEcB1wLxSb2Zml5jZGjNb89prr1W8sSIiIo2qZsl9ZnYPMDHk0JUx3+JzwM/d/SUzK3qiu98E3ASZOf4k7RQRKaXRSvSK5KpZ4Hf3s6KOmdkrZtbt7lvMrBt4NeS0k4BTzexzwFhglJltd/di+QAiIhXVqCV6RbLSMtS/DJgbPJ4LLC08wd0vdPfJ7j4FuAK4VUFfRGpNJXql0aUl8C8CzjazZ4Gzg58xs14z+25dWyYikkMleqXRpaKAj7v/ETgz5Pk1wMUhz98M3Fz1homIFJjU1cHmkCCvEr3SKNLS4xcRaQgq0SuNLhU9fhGRRqESvdLoFPhFRBKaM6NHgV4algK/iEgLU02C1qPALyLSolSToDUp8IuINKE4PfliNQkU+JuXAr+ISJOJ25NXTYLWpOV8IiJNJm51wajaA6pJ0NwU+EVEmkzcnrxqErQmBX4RkSYTtyc/Z0YPC887lp6uDgzo6epg4XnHan6/yWmOX0SkycybOS1vjh+ie/KqSdB6FPhFRJqMqgtKMQr8IiJNSD15iaI5fhERkRaiwC8iItJCFPhFRERaiAK/iIhIC1Fyn4iI5NGOfc1NgV9ERAZpx77mp6F+EREZFLfOvzQuBX4RERmkHfuanwK/iIgM0o59zU+BX0REBmnHvuan5D4RERmkOv/NT4FfRETyqM5/c9NQv4iISAtR4BcREWkhqQj8ZjbOzH5pZs8G3w+IOG+vmT0WfC2rdTtFREQaXSoCPzAfWOnuU4GVwc9hBtz9uODr3No1T0REpDmkJfDPBm4JHt8CzKljW0RERJpWWgL/BHffAhB8Hx9x3n5mtsbMHjazyJsDM7skOG/Na6+9Vo32ioiINKSaLeczs3uAiSGHrkzwNpPdvc/MDgdWmdkT7v584UnufhNwE0Bvb6+X1WAREZEmVLPA7+5nRR0zs1fMrNvdt5hZN/BqxHv0Bd83mtl9wAxgSOAXERGRcGkZ6l8GzA0ezwWWFp5gZgeY2ejg8YHAKcDTNWuhiIhIE0hL4F8EnG1mzwJnBz9jZr1m9t3gnOnAGjN7HLgXWOTuCvwiIiIJpKJkr7v/ETgz5Pk1wMXB44eAY2vcNBERkaaSlh6/iIiI1IC5N3fSu5m9Bvyh3u2ogAOB/6h3Iyqo2a4Hmu+amu16oPmuqdmuB5rvmup1PYe6+0FhB5o+8DcLM1vj7r31bkelNNv1QPNdU7NdDzTfNTXb9UDzXVMar0dD/SIiIi1EgV9ERKSFKPA3jpvq3YAKa7brgea7pma7Hmi+a2q264Hmu6bUXY/m+EVERFqIevwiIiItRIE/Zczsw2a2wcyeM7P5IcdHm9ntwfFHzGxK7VsZX4zr+aKZPW1m68xspZkdWo92JlHqmnLOu8DM3MxSldFbKM71mNnHg7+np8zsh7VuY1Ix/t1NNrN7zWxt8G/vo/VoZ1xm9j0ze9XMnow4bmb2z8H1rjOz42vdxiRiXM+FwXWsM7OHzOzdtW5jUqWuKee895rZXjO7oFZtG8Ld9ZWSL6CNzKZDhwOjgMeBdxac8zng28HjvwBur3e7h3k9pwNjgsefTfP1xL2m4Lw/Ax4AHgZ6693uYf4dTQXWAgcEP4+vd7srcE03AZ8NHr8TeKHe7S5xTacBxwNPRhz/KPB/AANOBB6pd5uHeT0n5/x7+0jaryfONQXntAGrgJ8DF9Srrerxp8v7gOfcfaO77wJ+DMwuOGc2cEvw+A7gTDOzGrYxiZLX4+73uvuO4MeHgYNr3Mak4vwdAVwLfA14s5aNK0Oc6/mvwL+6++sA7h66e2aKxLkmB94WPO4E+mrYvsTc/QFga5FTZgO3esbDQFew02kqlboed38o+++Nxvi9EOfvCOBvgTuJ2IG2VhT406UHeCnn503Bc6HnuPseoB94e01al1yc68n1aTK9ljQreU1mNgM4xN1/VsuGlSnO39GRwJFm9qCZPWxmH65Z68oT55quBj5pZpvI9L7+tjZNq5qk/9caSSP8XijJzHqAjwHfrndbUrFJjwwK67kXLruIc05axG6rmX0S6AU+UNUWDV/RazKzEcD1wEW1atAwxfk7GklmuP+DZHpevzKzY9x9W5XbVq441/SXwM3ufp2ZnQR8P7imfdVvXlU00u+F2MzsdDKB//31bksF3AB82d331nuQVoE/XTYBh+T8fDBDhyCz52wys5FkhilLDS/VS5zrwczOAq4EPuDuO2vUtnKVuqY/A44B7gv+c08ElpnZuZ7ZbTJt4v6be9jddwO/N7MNZG4E/r02TUwszjV9GvgwgLuvNrP9yNRUT/s0RpRY/9caiZm9C/gu8BHP7ODa6HqBHwe/Fw4EPmpme9x9Sa0boqH+dPl3YKqZHWZmo8gk7y0rOGcZMDd4fAGwyoOskRQqeT3BsPh3gHMbYO4YSlyTu/e7+4HuPsXdp5CZn0xr0Id4/+aWkEnCxMwOJDP0v7GmrUwmzjW9SLAVuJlNB/YDXqtpKytrGfCpILv/RKDf3bfUu1HlMrPJwF3Af3b3Z+rdnkpw98Nyfi/cAXyuHkEf1ONPFXffY2aXAivIZH9+z92fMrNrgDXuvgz4n2SGJZ8j09P/i/q1uLiY17MYGAv8JLgTftHdz61bo0uIeU0NI+b1rAA+ZGZPA3uBeWnugcW8pr8D/s3MLiczJH5Rim+gMbMfkZlqOTDIS/gHoB3A3b9NJk/ho8BzwA7gv9SnpfHEuJ6/J5O79K3g98IeT9lGN4ViXFNqqHKfiIhIC9FQv4iISAtR4BcREWkhCvwiIiItRIFfRESkhSjwi4iItBAFfhGpKDN70syursD7vGBmV1SgSSKSQ+v4RaSugpuEC9z9mIJD7wXeqH2LRJqbAr+IpJK7N3IlPZHU0lC/SIMys/vM7Ntm9k0zez34WhxsFISZnWdm68xswMy2mtn9ZjYh5/WzzOw3Zvammf3ezP4xKHGbPT5kqD34zBtzfh5vZkuDz/iDmf11SDsnm9lPzez/BV93mdnBwbGLyFQ4O9rMPPi6KOzzg2OfDT5vh5k9Y2anm9nBZrbCzN4ws8fM7PiCzz85uPYdZrbZzP6Hmb0NkRalwC/S2C4k8//4JOC/AZcAl5nZRDL70N8CTAdOA76ffZGZzQRuA24Ejgb+mszeD/+U8PNvBt4BnAXMAT4FTMn5HCNT638CcAaZmv+TgCXBsduB64ANQHfwdXuRz7squK53A2uAH5EpY/0tYAaZjWluzvn8Y4FfkKll/27gPOA44HsJr1OkaWioX6SxbQE+H9SZ/52ZHQl8EbiPTJ3wO9z9D8G5T+a87kpgsbv/r+Dn583sy8APzGxenLr1wWd9BHi/uz8YPDeX/A18ziITcI9w9xeCc/6KTE35M939HjPbTqYW+8sxrvdWd/9R8D7/RGZ73RXuvjR47mvAvWZ2oLv/BzAPuN3dr8tp92eBtWY2vkE2hhKpKPX4RRrbwwVBejXQAzwP3AM8aWZ3BkPkB+Wc9x7gSjPbnv0CfgjsT2Yr4TimA/uAR7NPBDcZfQXn9GWDfnDOxuCcd8b8nFzrch6/Enx/IuS58cH39wCfLLjOB4NjR5Tx+SINTz1+kebkwIeAE4PvnwYWmtkH3P1xMjf9/x34Schrs0l1+wArONae87jwWBgL2hLVxqR2h7w+7LkROd+/C1wf8l6by/h8kYanwC/S2E4wM8vp9Z9Ipof9p+Dn1cDqYEvap4BPAI8DvwWOcvfnirz3a2Tm3AEws/2Ao4C1wVPryQTW9wIPBedMJjOHn/U00GNmU3KG+g8Pznk6OGcXme1zq+G3wNElrlOkpWioX6SxTQJuMLNpZnYBmTnt683sRDO7yszeGwTjc4FDeCvYXgP8lZldY2bHmNlRZnZBMEeetQq40Mw+aGZHk0mIG+zxu/sG4P8C3zGzk8zsODKJdQM573EPmRuN28zsPWbWSyap8LfB+wO8ABxqZseb2YFmNrpyfzx8FXhfsPphhpm9w8z+k5l9p4KfIdJQFPhFGtttZHrLjwD/RibD/XqgHzgF+BnwLJnM+Wvd/QcA7r4COIdMlv2jwdd84MWc915IJjgvJZMZ/2syATvXRcDvg/PuJpMn8EL2YDASMYfM6MF9wL3Ay8CcnFGKO4GfAyuD8/6yzD+LIdx9HZkVDVOA+8nchCzkrVwAkZZjMZJ3RSSFzOw+4El3v7TebRGRxqEev4iISAtR4BcREWkhGuoXERFpIerxi4iItBAFfhERkRaiwC8iItJCFPhFRERaiAK/iIhIC1HgFxERaSH/H8/G47r7Ly7YAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8, 6))\n",
    "colors = plt.get_cmap(\"tab10\").colors[::-1]\n",
    "labels = df.index.unique()\n",
    "\n",
    "X = gplvm.X_loc.detach().numpy()\n",
    "for i, label in enumerate(labels):\n",
    "    X_i = X[df.index == label]\n",
    "    plt.scatter(X_i[:, 0], X_i[:, 1], c=[colors[i]], label=label)\n",
    "\n",
    "plt.legend()\n",
    "plt.xlabel(\"pseudotime\", fontsize=14)\n",
    "plt.ylabel(\"branching\", fontsize=14)\n",
    "plt.title(\"GPLVM on Single-Cell qPCR data\", fontsize=16)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the first dimension of the latent $X$ for each cell (horizontal axis) corresponds well with the observed capture time (colors). On the other hand, the 32 TE cell and 64 TE cell are clustered near each other. And the fact that ICM cells differentiate into PE and EPI can also be observed from the figure!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Remarks\n",
    "\n",
    "+ The sparse version scales well (linearly) with the number of data points. So the GPLVM can be used with large datasets. Indeed in [2] the authors have applied GPLVM to a dataset with 68k peripheral blood mononuclear cells.\n",
    "\n",
    "+ Much of the power of Gaussian Processes lies in the function prior defined by the kernel. We recommend users try out different combinations of kernels for different types of datasets! For example, if the data contains periodicities, it might make sense to use a [Periodic kernel](http://docs.pyro.ai/en/dev/contrib.gp.html#periodic). Other kernels can also be found in the [Pyro GP docs](http://docs.pyro.ai/en/dev/contrib.gp.html#module-pyro.contrib.gp.kernels)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References\n",
    "\n",
    "[1] `Resolution of Cell Fate Decisions Revealed by Single-Cell Gene Expression Analysis from Zygote to Blastocyst`,<br />&nbsp;&nbsp;&nbsp;&nbsp;\n",
    "Guoji Guo, Mikael Huss, Guo Qing Tong, Chaoyang Wang, Li Li Sun, Neil D. Clarke, Paul Robson\n",
    "\n",
    "[2] `GrandPrix: Scaling up the Bayesian GPLVM for single-cell data`,<br />&nbsp;&nbsp;&nbsp;&nbsp;\n",
    "Sumon Ahmed, Magnus Rattray, Alexis Boukouvalas\n",
    "\n",
    "[3] `Bayesian Gaussian Process Latent Variable Model`,<br />&nbsp;&nbsp;&nbsp;&nbsp;\n",
    "Michalis K. Titsias, Neil D. Lawrence\n",
    "\n",
    "[4] `A novel approach for resolving differences in single-cell gene expression patterns from zygote to blastocyst`,<br />&nbsp;&nbsp;&nbsp;&nbsp;\n",
    "Florian Buettner, Fabian J. Theis"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
